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The  aspects  of  accurate  determination  of  Earth  satellite  orbits  by  the 
Minimum  Variance  Method  are  presented.  In  addition,  techniques  for  the 
determination  of  the  associated  physical  constants,  such  as  the  coefficients 
in  the  Earth's  gravitational  potential,  exospheric  temperature,  etc. .  are 
developed.  A  method  for  determination  of  the  state  transition  matrix  is 
presented.  Also  include  J  are  a  review  of  the  time  systems  employed  in 
satellite  orbit  determination  and  a  short  discussion  of  the  types  of  observa¬ 
tions. 

The  mathematical  model  of  the  dynamical  system  includes  nine  zonal  har¬ 
monics  and  up  to  the  fourth  order  tesseral  harmonics  of  the  Earth's  gravi¬ 
tational  potential.  Atmospheric  drag  effects  are  included  on  the  assumption 
that  the  atmosphere  rotates  v/ith  the  angular  velocity  of  the  Earth.  First 
order  solar  and  lunar  gravitational  attractions  and  solar  radiation  pressure 
arc  also  treated.  The  satellite  orbits  are  integrated  in  a  reference  system 
which  considers  the  precession  and  nutation  of  the  Earth.  Rectangular  co¬ 
ordinate  systems  are  used  throughout  the  development. 
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FOREWORD 

This  report  presents  the  analytical  foundations  for  several  computer  pro¬ 
grams  now  under  development.  The  research  is  being  performed  for  the  Data 
Analysis  Branch  (CRMXA),  Technical  Services  Division  at  AFCRL,  USAF, 

L.  G.  Hanscom  Field,  Bedford,  Massachusetts.  The  digital  computer  pro¬ 
gramming  in  this  research  effort  is  being  done  by  Bruce  Clemenz  and  Jacques 
Fein.  The  contractor's  report  number  is  ER  13950. 
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L  W/Y.'  FICTION 

The  orbit  of  a  satellite,  in  a  given  refc*-  ;;e  system  is  determined  by  six 
parameters,  called  the  oruital  elements.  These  elements  ~an  be  determined 
from  a  sufficient  number  of  other  parameters  obtained  a.’  observations.  How¬ 
ever.  the  problem  is  usually  either  under-  or  over-determined.  In  adoition, 
the  observations  contain  errors.  For  these  reasons,  normal  algebraic  .  lethods 
cannot  be  employed,  and  the  problem  enters  the  realm  of  statistics  and  prob¬ 
ability.  Satellite  observations  are  used  not  only  to  determine  the  orbit,  but 
also  to  estimate  the  various  physical  and  other  constants  of  the  dyn:  nical  and 
the  observation  systems.  Precise  orbit  determination  and  analysis  uy  a  statis¬ 
tical  filtering  technique  known  as  the  Minimum  Variance  Method  is  the  subject 
of  this  report. 

The  motion  of  a  satellite  in  an  inverse -square  central  force  fielu  is  repre¬ 
sented  by  conic  sections.  This  motion  would  be  realized  if  the  sole  force  act¬ 
ing  cn  the  satellite  were  due  to  a  point  mass  or  a  homogeneous  sphere.  In 
nature,  however,  the  orbit  of  a  satellite  is  perturbed  by  a  variety  of  forces, 
e.g. ,  harmonics  in  the  gravitational  potential,  atmospheric  drag,  attractions 
of  other  celestial  bodies,  etc.  Consequently,  the  dynamical  system  is  much 
more  complicated. 

As  mentioned  before,  the  orbit  of  a  satellite  is  determined  by  the  six  orbi¬ 
tal  elements.  These  elements  can  represent  an  orbit  in  a  simple  inverse- 
square  central  force  field  as  well  as  in  the  actual,  complex  dynamical  system. 
Any  set  of  independent  parameters  which  uniquely  describe  an  orbit  can  be  used 
as  orbital  elements.  For  example,  the  six  classical  elements  are:  a,  the  semi¬ 
major  axis;  e,  the  eccentricity;  ft,  the  right  ascension  of  the  ascending  node; 

<a,  the  argument  of  pericenter;  i,  the  inclination;  and  M,  the  mean  anomaly. 

The  orbits  in  the  present  report  are  considered  in  a  rectangular  coordinate 
system.  Consequently,  the  orbital  elements  in  this  system  are  most  conven¬ 
iently  represented  by  the  six  components  of  the  position  and  velocity  vectors, 
which  are  the  standard  elements  in  this  analysis.  The  certain  numerical  dif¬ 
ficulties  experienced  with  these  elements  in  the  Least  Squares  application  are 
not  encountered  employing  the  Minimum  Variance  technique. 

The  standard  reference  epoch  for  the  celestial  equations  in  this  report  is 
1950  January  1,0  UT,  except  as  noted. 

The  present  report  deals  specifically  with  the  orbit  determination  of  close 
Earth  satellites,  including  parabolic  and  hyperbolic  orbits.  The  same  princi¬ 
ple,  however,  applies  to  satellites  of  other  celestial  bodies. 
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II.  TIME  SYSTEMS 


A.  MEASUREMENT  OF  TIME 

The  accurate  measurement  of  time  has  been  the  concern  of  a  lather  spec¬ 
ialized  field.  However,  with  the  advent  of  the  artificial  satellite,  a  new  area 
has  been  introduced  which  is  directly  concerned  with  the  px'ecise  aspects  of 
time  and  its  measurement.  In  addition,  there  have  been  many  developments 
in  the  time  systems  during  the  last  few  years.  Therefore,  a  review  of  this 
subject  is  necessary.  A  thorough  treatment  can  be  found  in  Refs.  1,  2  and  3. 
The  present  discussion  will  deal  only  with  the  existing  time  systems  and  their 
application  to  satellite  orbit  determination  and  analysis. 

In  the  measurement  of  time,  two  important  factors  are  involved.  One  is 
the  epoch  or  reference  from  which  time  is  measured.  The  otl.°r  is  the  rate  at 
which  time  is  measured.  A  system  used  to  measure  time,  therefore,  must 
be  some  observable  pnysical  phenomenon  that  has  a  well-defined  epoch  and 
whose  rate  of  change  is  as  invariant  as  possible.  There  are  seven  fundamental 
time  systems  in  existence:  ephemeris  time  (ET),  atomic  time  (A.  1),  true 
sidereal  time,  mean  sidereal  time,  UTO,  UT1  and  UT2.  The  last  three  are 
subdivisions  of  universal  time.  In  addition,  there  are  time  systems  emitted 
as  signals  or  several  radio  frequencies  throughout  the  world:  WWV,  NBA 
(United  States) ,  GBR  (United  Kingdom) ,  etc. 


B.  EPIIEMERIS  TIME 

Ephemeris  time  is  the  uniform  time  of  dynamical  astronomy.  Theoretically, 
it  is  the  independent  time  argument  of  the  ephemerides  of  the  Sun,  Moon  and 
the  planets.  In  practice,  it  is  determined  from  the  orbital  motion  of  the  Moon 
about  the  Earth.  This  involves  the  solution  of  the  equations  of  motion  of  the 
Moon.  Thus,  the  uniformity  of  the  time  argument  will  be  dependent  on  the  ac¬ 
curacy  of  the  representation  of  the  system  in  which  the  Moon  moves.  Since  it 
also  involves  observations  of  the  Moon,  additional  errors  are  introduced.  The 
errors  are  not  significant,  and  the  accuracy  of  the  determined  ET  is  believed 
to  be  within  a  couple  of  seconds  in  a  century. 

Ephemeris  time  would  be  the  proper  time  argument  in  the  equations  of  mo¬ 
tion  of  an  Earth  satellite.  However,  its  determination  requires  several  years 
of  observations  and  it  is  not  practical  for  observing  events  of  relatively  short 
time  intervals  apart  without  an  intermediary. 


C.  ATOMIC  TIME 

In  1955,  a  precise,  cesium-133  atomic  resonator  was  introduced  in  Great 
Britain.  Since  then,  nine  laboratories  in  five  countries  have  been  operating 
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cesium  resonators  in  a  coordinated  effort.  It  is  estimated  that  the  accuracy  of 
the  cesium  resonators  is  within  0.001  second  in  three  years.  Because  of  the 
exceptional  accuracy  of  this  system  (designated  A.  1),  the  12th  General  Con¬ 
ference  on  Weights  and  Measures,  in  December  1964,  changed  the  basis  of  the 
definition  of  the  second  from  ephemeris  to  atomic  time  scale.  Thus,  a  uniform 
time  system  has  been  made  available  which  provides  exceptional  stability  and 
convenience. 

The  A.  1  system  is  defined  as  follows  (Ref.  2): 

(1)  A  clock  which  keeps  A.  1  time  advances  one  second  in  the  interval 
required  for  9, 192,631,770  oscil.ations  of  cesium  at  zero  field. 

(2)  At  0h  oV  UT2  on  1  January  1958,  the  value  of  A.  1  was  0h  0™ 0s. 

Because  the  atomic  time  system  is  uniform  for  all  practical  purposes  and 
readily  available,  it  is  the  most  suitable  for  use  as  the  independent  time  argu¬ 
ment  in  the  satellite  equations  of  motion. 

The  rapid  progress  in  atomic  timekeeping  is  indicated  by  the  announcement, 
even  as  this  report  was  being  written,  that  the  U.  S.  Naval  Research  Labora¬ 
tory  has  installed  twin  atomic  hydrogen  masers  to  continuously  reset  the  master 
clock  at  the  U.  S.  Naval  Observatory.  The  hydrogen  masers  will  keep  the  clock 
accurate  to  within  one  second  in  300,000  years. 


D.  SIDEREAL  TIME 

Sidereal  time  is  directly  related  to  the  rotation  of  the  Earth.  It  is  defined 
as  the  hour  angle  of  the  vernal  equinox.  Thus,  except  for  small  motions  of  the 
equinox  itself,  sidereal  time  is  a  direct  measure  of  the  diurnal  rotation  of  the 
Earth.  Sidereal  t;me  measured  with  respect  to  the  true  equinox  is  true  or 
apparent  sidereal  time.  If  measured  with  respect  to  the  mean  equinox  of  date, 
it  is  called  mean  sidereal  time.  If  the  meridian  is  that  of  Greenwich,  it  is 
called  Greenwich  mean  sidereal  time.  For  any  other  meridian,  it  is  the  local 
sidereal  time.  Because  sidereal  time  reflects  the  variable  rotation  of  Earth, 
there  is  no  direct  relationship  between  sidereal  time  and  ephemeris  time.  In 
orbit  determination,  the  computation  of  sidereal  time  is  required  in  order  to 
determine  the  position  of  the  observing  station. 


E.  UNIVERSAL  TIME 

Universal  time  is  based  on  the  diurnal  motion  of  the  Sun  and  is  the 
basis  of  all  civil  timekeeping.  It  involves  both  the  rotation  of  Earth  and  the 
motion  of  Earth  in  its  orbit  about  the  Sun.  Universal  time  and  sidereal  time 
are  equivalent  systems  and  are  directly  related  to  each  other  by  means  of  a 
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numerical  formula.  Universal  time  is  defined  as  the  Greenwich  hour  angle  of 
a  point  on  the  equator  whose  right  ascension,  measured  from  the  mean  equinox 
of  date,  is 
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Time,  Tu,  is  given  in  Julian  centuries  of  36,525  days  of  universal  time  elapsed 

since  the  epoch  1900  January  0,  12h  UT.  The  practical  determination  of  univer¬ 
sal  time  is  made  through  the  intermediary  of  sidereal  time  by  observing  the 
diurnal  motions  of  the  stars.  The  relation  between  Greenwich  mean  sidereal 
time  and  UT  is  given  by  the  equation: 

Greenwich  mean  sidereal  time  -  UT  -*  Ru  +  12h  (2) 

The  sidereal  times  are  computed  in  advance  for  successive  dates  of  0h  UT 
from  the  above  equations  and  published  in  the  American  Ephemeris  and  Nautical 
Almanac.  Local  mean  sidereal  time  at  any  particular  instant  is  obtained  from 
observations  of  the  transit  of  stars  of  known  positions.  Greenwich  mean  side¬ 
real  time  is  obtained  by  adding  the  lor.gitudc  west  of  Greenwich.  Then  the  cor¬ 
responding  universal  time  is  obtained  by  taking  the  difference  in  sidereal  times 

at  the  instant  of  the  observation  and  the  computed  value  at  0h  UT  and  converting 
it  to  UT  by  the  known  relationship.  The  universal  time  thus  obtained  is  desig¬ 
nated  UTO.  Due  to  the  movement  of  the  Earth's  axis  of  rotation,  known  as 
polar  motion,  and  variation  in  the  rate  of  rotation  of  the  Earth,  UTO  contains 
irregularities.  Although  the  polar  motion  is  very  small,  it  affects  the  time 
measurements,  which  can  be  now  performed  with  great  accuracy.  UTO  cor¬ 
rected  for  the  polar  motion  to  a  mean  Greenwich  meridian  is  designated  UTl, 


The  variations  in  the  rate  of  rotation  of  the  Earth  are  due  to  many  causes. 
Some  of  them  are  negligible;  some  of  them  are  quite  large  but  unpredictable. 
The  seasonal  variation  is  a  periodic  variation  of  a  maximum  amplitude  of  ap¬ 
proximately  0.03  second.  It  is  a  quite  stable  variation  and  can  be  predicted 
with  good  accuracy.  UTl  corrected  for  the  seasonal  variations  is  designated 
UT2.  Corrections  for  the  years  1956  to  1962  have  been  computed  from  the 
formula  (Ref.  1) : 

AT  =  +  0s.  022  sin  2m  -  0*.  017  cos  2ttt 

s  s  <3> 

-  0.007  sin  4ttt  +  0. 006  cos  4ttt 

Since  1962,  the  following  formula  has  been  used  (Ref.  4): 

AT  -  +  0S.  022  sin  2irr  -  0*012  cos  2m 

s  s  (4) 

-  0.006  sin  4ttt  +  0.007  cos  4ttt 


where  t  is  the  fraction  of  the  year  and  is  zero  on  January  1. 


It  is  important  to  note  that  universal  time  is  not  a  uniform  time  and  therefore 
cannot  be  properly  used  as  the  independent  time  argument  in  the  equations  of 
motion.  However,  UT1  is  of  particular  significance  in  accurate  satellite  orbit 
determination.  Since  UT1  is  obtained  by  observing  the  rotational  position  of 
Earth  with  respect  to  stars,  the  reverse  process  is  applied  to  determine  the 
position  of  Earth  from  UTl. 

As  noted  before,  only  the  ratio  of  sidereal  time  and  universal  time  can  be 
expressed  by  a  numerical  formula.  There  are  no  rigorous  analytical  relation¬ 
ships  for  the  other  systems.  The  difference  between  ephemeris  time  and 
universal  time  presently  is  about  35  seconds.  The  difference  between  A.l  and 
UT2  was  2.3385  seconds  on  January  0,  1963,  increasing  by  about  0.  5  second 
per  year. 


F.  RADIO  TIME  SIGNALS 

Neither  of  the  time  systems  discussed  above  is  directly  available  to  the  user. 
Instead,  time  signals  are  emitted  by  special  radio  stations  (WWV,  WWVH  and 
NBA  in  the  United  States)  in  accordance  with  international  agreements.  Prior 
to  1959,  station  WWV  emitted  time  signals  at  a  constant  frequency  while  making 
phase  adjustments  of  0.  02  second  when  necessary  to  keep  the  signals  close  to 
UT2.  The  signals  now  are  emitted  with  a  frequency  maintained  constant  each 
year  but  offset  with  reference  to  the  atomic  time  standard.  The  time  pulses 
are  kept  within  approximately  0. 1  second  of  UT2.  In  addition,  phase  adjust¬ 
ments  of  the  pulses  can  be  made  if  necessary.  There  is  no  analytical  relation¬ 
ship  between  WWV  and  UT2.  For  this  reason,  the  differences  are  given  in 
periodical  Time  Service  Bulletins  issued  by  the  U.  S.  Naval  Observatory. 

The  methods  used  in  satellite  observation  and  timing  are  continuously  im¬ 
proving.  It  is  estimated  (Ref.  5)  that  the  position  of  a  satellite  can  be  meas¬ 
ured  with  an  accuracy  of  20  meters.  For  practical  reasons,  the  timing  of  a 
satellite  observation  is  done  by  the  clock  of  a  station  which  observes  the  satel¬ 
lite.  Although  the  received  WWV  signals  will  be  in  error  with  respect  to  the 
emitted  signals  because  of  uncertainty  in  propagation,  it  is  estimated  (Ref.  5) 
that  a  worldwide  tracing  system  can  be  synchronized  to  an  accuracy  of  about 
0.001  second.  To  fully  utilize  the  accuracy  available  in  precision  orbit  determina¬ 
tion,  the  proper  corrections  should  be  applied  to  the  observation  times  to  arrive 
at  the  correct  time  systems  required  in  the  analysis.  These  corrections,  how¬ 
ever,  are  not  immediately  available  and  in  such  cases  they  must  be  either 
extrapolated  or  the  time  recorded  by  the  station  clock  used  as  an  approxima¬ 
tion.  During  the  period  in  which  there  are  no  changes  either  in  frequency 
or  phase  of  the  WWV  signals,  the  recorded  time  will  be  essentially  uniform, 
if  the  station  clocks  are  well  synchronized  with  the  WWV.  Therefore,  this 
time  could  be  used  as  the  independent  time  argument  in  the  satellite  equations 
of  motion.  An  error  is  introduced  by  substituting  this  time  for  UTl  to  compute 
the  station  position.  For  example,  the  position  of  a  station  at  30°  latitude  can 
be  in  error  by  as  much  as  50  meters.  The  numerical  values  for  the  masses  of 


Earth  and  other  celestial  bodies  are  not  sufficiently  well  known  at  the  present 
time  to  be  affected  by  the  small  differences  in  the  time  systems,  but  this  may 
change  in  the  future. 

The  following  procedure  can  be  used  to  obtain  systems  A.  1  and  UTl  required 
for  accurate  orbit  analysis: 

(1)  Time  recorded  by  the  station  clock  corrected  to  obtain  WWV  emitted 
signal 

(2)  WWV  emitted  signal  corrected  to  obtain  A.  1 

(3)  WWV  or  A.  1  corrected  to  obtain  UT2 

(4)  UT2  corrected  to  obtain  UTl. 

Correction  (1)  includes  eorreetions,  if  any,  to  the  recorded  time  to  obtain  WWV 
received.  This  time  is  then  corrected  for  propagation  effects  to  obtain  WWV 
emitted.  Corrections  (2),  (3)  and  (4)  are  published  in  the  Time  Service  Bul¬ 
letins.  Correction  (4)  can  be  applied  by  means  of  Eqs.  (3)  and  (4). 


III.  SPACE  REFERENCE  SYSTEMS 


A.  SYSTEM  REQUIREMENTS 

There  are  several  basic  reference  systems  used  in  orbit  determination 
The  satellite  motion  itself  must  be  ultimately  considered  in  an  inertial  sys¬ 
tem.  The  position  of  the  observer  must  be  referenced  to  a  terrestrial  sys¬ 
tem,  and  the  observations  are  obtained  in  either  a  geodetic  or  a  celestial 
system.  In  the  final  analysis,  the  relationship  between  these  systems  must 
be  introduced. 


B.  BASIC  REFERENCE  SYSTEM 

The  satellite  orbit  could  be  referenced  to  a  fixed  geocentric  celestial  sys¬ 
tem.  However,  such  a  system  would  not  be  very  convenient  for  determining 
the  forces  acting  on  the  satellite  because  of  the  precession  and  nutation  of  the 
Earth.  The  expressions  for  the  forces  could  be  much  simplified  if  the  motion 
is  considered  in  a  moving  axis  system  defined  by  the  true  equator  and  equinox. 
In  such  a  system,  however,  a  supplementary  or  Coriolis  acceleration  is  intro¬ 
duced.  Since  the  rate  of  precession  and  nutation  is  comparatively  small,  we 
will  utilize  a  system  which  is  considered  inertial  for  a  short  period  of  time 
and  coincides  with  a  mean  position  of  the  true  equator  and  equinox  during  this 
interval.  The  satellite  orbit  thus  is  considered  in  a  system  which  moves  step¬ 
wise  with  the  rate  of  precession  and  nutation.  The  length  of  the  interval  can 
be  made  as  short  as  is  necessary  and  practical.  The  only  errors  introduced 
will  be  due  to  the  small  variation  of  the  gravitational  field  caused  by  preces¬ 
sion  and  nutation  during  this  interval.  However,  for  all  practical  purposes, 
these  errors  will  be  negligible.  In  addition  to  giving  minimum  errors,  this 
system  is  ideally  suited  to  the  point-to-point  technique  of  the  Minimum  Vari¬ 
ance  Method  and  the  integration  method  employed. 


C  .  TERRESTRIAL  REFERENCE  SYSTEM 

A  terrestrial  reference  system  will  be  defined  as  a  rectangular  right-hand 
system  with  the  origin  at  the  Earth's  center  of  gravity  and  the  z-axis  directed 
toward  the  mean  north  pole  as  defined  by  the  International  Latitude  Service. 

The  x-z  plane  will  coincide  with  the  mean  meridian  of  Greenwich.  This  defini¬ 
tion  is  purely  theoretical,  since  the  Earth's  center  of  gravity  is  not  precisely 
known.  For  this  reason,  the  positions  given  in  this  system  will  be  in  error  and 
thus  affect  the  satellite  observations.  A  method  to  improve  observing  station 
positions  from  orbit  analysis  by  the  Minimum  Variance  Method  will  be  presented 
later  in  the  report. 
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Generally,  the  coordinates  of  an  observing  station  are  given  as  polar 
coordinates  in  a  gcodetie  system.  The  coordinates  are  longitude,  latitude 
and  altitude.  Longitude  will  be  considered  positive  west  for  the  station  lo¬ 
cations.  Latitude  is  usually  given  as  geodetie  or  geographic  latitude.  The 
two  differ  due  to  local  gravity  anomalies.  Assuming  that  the  given  latitude 
is  geodetic,  we  must  obtain  the  geocentric  latitude  as  a  first  step  toward 
conversion  to  rectangular  coordinates.  It  is  important  to  note  that  the  given 
polar  coordinates  are  always  associated  with  a  speeifie  ellipsoid,  defined  by 
the  mean  equatorial  radius  of  the  Earth,  R  ,  and  flattening,  f.  The  values  for 

the  adopted  International  Ellipsoid  of  Reference  are  R,.  -  6,378,388  meters,  and 

f  1/297.  More  recent  determinations  have  given  R^  =  6,378, 156  meters  and 

f  --  1/29S.3. 

The  transformation  from  geodetic,  T ,  to  geocentric  latitude,  ean  be 
obtained  by  considering  the  geometry. 


Fig.  1.  Station  Coordinates 


Designating  f  -  (1  -  f)2  we  ean  write: 


tan  - 


(p‘  +  h)  sin  t 
h  eos  ■+  eos  <{> 


(5) 


where  p( ,  p'  and  h  are  expressed  in  units  of  R^.  Introducing  an  auxiliary 
function,  C,  defined  by: 


p(  cos  <t>,  C  eos  ? 
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/  2  \1/- 
cos  t  [1  +  fF  tan 

and  since  p  =  f,.C,  we  can  express  4>  as  function  of 

1  t, 

h  +  f  C 

L. 

tan<£  =  ITTc  tan 

Then  the  geocentric  radius  of  the  observing  station  is: 


P  = 


Re  (h  +  C) 


cos  f 

COS  c{> 


(C) 


(7) 


The  station  coordinates  in  the  Earth  fixed  terrestrial  system  are 


Xg'  =  p  cos  <f>  cos  (-X) 

(8a) 

y”  =  p  cos  sin  (-X) 

(3b) 

Zg  =  p  sin  <)> 

(8c) 

where  X  is  the  station  longitude  (positive  west). 

These  coordinates  are  in  the  terrestrial  system  as  previously  defined. 

The  actual  Earth’s  axis  of  rotation,  however,  does  not  coincide  with  the  mean 
axis  as  defined  by  the  International  Latitude  Service.  It  is  moving  about  the 
mean  pole  in  what  is  known  as  the  polar  motion.  This  motion  has  the  effect 
of  a  small  variation  in  latitude  and  meridian  and  must  be  considered  in  ac¬ 
curate  calculations.  The  variations  are  regularly  published  by  the  Interna¬ 
tional  Latitude  Service.  The  transformation  of  the  station  coordinates  from 
the  mean  to  the  instantaneous  system  is  accomplished  by  a  simple  transforma¬ 
tion: 


where  x  and  y  are  the  angular  coordinates  of  the  instantaneous  pole,  in  radians. 


D.  STATION  POSITION  IN  THE  BASIC  SYSTEM 

So  far  the  obtained  position  of  the  observing  station  is  in  a  rotating  Earth 
fixed-axes  system.  To  relate  it  to  the  previously  defined  basic  reference 


11 


system,  which  lor  all  practical  purposes  is  a  true  sidereal  system,  one  more 
transformation  must  be  performed. 

The  Earth  fixed  system  and  the  basic  system  have  the  same  z-axis  (axis  of 
rotation).  The  position  of  the  observing  station  in  the  basic  system,  therefore, 
can  be  obtained  by  utilizing  the  true  sidereal  time,  which  is  a  function  of  the 
universal  time  (UT1).  This  will  give  the  Greenwich  hour  angle  of  the  vernal 
equinox  as  computed  from  th  i  dlowing  equation: 

X  =  1.  74GG47719  +  0.30038809863056  d  +  0.5064  x  10'1"1  d2+  AA  (10) 

where  Af;Rand  AAare  in  radians  and  d  in  Julian  days  from  the  epoch  1950 

January  1,  0h  UT.  As  discussed  in  Chapter  II,  d  must  be  expressed  in  UT1. 
The  above  equation  is  based  on  Newcomb's  expression  and  corresponds  to 
values  for  sidereal  times  published  in  the  American  Ephemeris  and  Nautical 
Almanac  from  1960  on.  The  terms  due  to  nutation  are  expressed  by  A  A. 

AA  -  -  0.  76700  x  10~4  sin  (0.  211408  -  0.00092422  d) 

+  0.929  x  10'f>  sin  (0.422816  -  0.00184844  d) 

-  0.907  x  10'fl  sin  (2.247127  *  0.45994300  d)  (11) 

-  0.5662  x  10' 5  sin  (9.776679  4  0.03440558  d) 

0.560  x  10'f>  sin  (G.  248291  4  0.01720197  d) 


The  station  coordinates  in  the  basic  system  are  then: 


Xs 

o 
o 
c n 

V 

sinAc;R 

0 

1 

X 

*» 

y* 

= 

sin  AC;R 

cos  Af.R 

0 

1 

ys 

(12) 

Zs 

0 

0 

1 

i 

Z- 

a 

E.  PRECESSION  AND  NUTA  ITON 


The  basic  system,  as  defined  previously,  is  essentially  a  true  sidereal  sys¬ 
tem.  As  such,  it  follows  the  precessional  and  nutational  motions  of  the  Earth. 
The  coordinate  transformation  from  the  reference  mean  equinox  of  1950 

January  1,  0h  UT  to  the  mean  equinox  of  Hate  is  arrnmplished  by  the  following 
matrix : 


P 


X 


V 


X, 


Yv 

Y. 


Zv 

z . 


(13) 
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where-. 


Yn 

Z, 

Yz 


„t9  9  -20 

Xx  =  1  -  0.2226036  x  10  "  d  -  0.267  x  10  d 


-Xy  - 

-  0.61190636  x  10'h 

d  -  0.5067  x  10'14 

d2  + 

0.453  x 

io'! 

d 

-G 

- 14 

2 

0. 197  x 

-l 

1!'  ■ 

-x,  - 

-  0.26603997  x  10 

d  +  0. 1552  x  10 

d  + 

10 

d 

.v)  2  -20  'i 

Yv  =  1  -  0.1872158  x  10  d“  -  0.308  x  10  dJ 

Zv  =  -  0.  813972  x  10‘13  d2  -  0.  61  x  10'21  d3 

Zz  =  1  -  0.353878  x  1013  d2  +  0.41  x  10'21  d3 


and  d  are  Julian  days  since  1950  January  1,  0  UT.  To  account  for  the  nutation, 
a  transformation  from  the  mean  equinox  of  date  to  the  true  equinox  of  date  is 
accomplished  by  the  following,  matrix: 


1 

-Ap 

-Av 

Ap 

1 

-Ae 

Av 

Ae 

1 

(14) 


where  Ap,  Av  and  Ac  are  the  terms  due  to  nutation  in  right  ascension,  declina¬ 
tion  and  obliquity,  respectively.  They  can  be  computed  with  sufficient  accuracy 
from  the  following  expressions 

Ap  =  (-76.  700  sir.  ol  +  0.  929  sin  o2  -  0.  907  sin  -  5.  662  sin  o-4 
+  0.560  sin  or )  x  10  b 

Av  =  (-33.  009  sin  <*l  +  0.400  sin  a>2  -  0.390  sin  a3  -  2.  437  sin  r*4 
+  0. 241  sin  a5)  x  10‘6 

Ae  =  (+44.  654  cos  -  0.438  cos  a2  +  0.428  cos  a3  +  2.676  cos  a4)  x  10'6 
where 

=  0.211408  -  0.00092422  d 
u2  =  0.422816  -  0.00184844  d 
«3  =  2.247127  +  0.45994300  d 

=  .9.776679  +  0.03440558  d 
o5  =  6.248291  +  0.01720197  d 
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Only  the  most  significant  nutation  terms  arc  considered.  However,  they  will 
give  an  accuracy  of  better  than  four  meters  on  the  Earth's  surface. 

The  transformation  of  the  rectangular  coordinates  between  two  basic  or 
true  sidereal  systems  of  arbitrary  dates  can  be  accomplished  by  the  following 
matrix  multiplication: 

x,  -  N,  P,  P,r  N[  xt]  (15) 

where  _  „  - 


r  i 

• 

* 

X 

* 

>• 

or 

y 

* 

_z . 

t 

z 

and  T  means  the  transpose  of  a  matrix. 

In  a  continuous,  point-to-point  transformation,  only  three  multiplications  must 
be  performed  at  every  transformation,  since  the  transpose  of  the  previous  N2  P2 

T  T 

matrix  is  the  new  P(  N(  matrix. 

It  is  convenient  to  consider  the  precession  and  nutation  matrices  as  sums 
of  a  unit  matrix  and  a  matrix  with  small  elements. 

Designating 

P '  P  -  1 


and 


N'  N  -  1 

where  I  is  a  unit  matrix. 

NP  =  (N ’  -  I)  (P 1  +  I)  -  N'P’  +  P’  +  N'  +  I 

Matrices  N'P',  P'  and  N'  contain  only  small  elements.  Consequently,  fewer 
significant  figures  need  to  be  carried.  The  number  of  significant  figures  will 
increase  for  the  diagonal  elements  only  in  the  addition  to  the  unit  matrix.  In 
some  cases,  it  might  be  possible  to  neglect  the  second  order  term  N'P'  and 
thus  the  matrix  multiplication  could  be  avoided  completely.  Since  the  preces- 
sional  and  nutational  motion  is  of  the  order  of  approximately  0'.'3  per  day,  the 
transformation  need  not  be  performed  very  often. 
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F.  OBSERVATION  REFERENCE  SYSTEMS 


Two  types  of  observations  will  be  considered:  (I)  observations  referenced 
to  a  geodetic  system  and  (2)  observations  referenced  to  a  celestial  system. 

Both  types  are  necessarily  topocentr’c  systems.  Types  of  observations  refer¬ 
enced  to  a  geodetic  system  will,  generally,  include  elevation,  azimuth,  eleva¬ 
tion  and  azimuth  rate,  range,  range  rate  and  range  acceleration.  Obse na¬ 
tions  referenced  to  a  celestial  system  are  right  ascension  and  declination. 

To  perform  the  orbit  analysis,  the  observations  corresponding  to  the  estimated 
orbit  must  be  compared  with  the  actual  observations  in  a  common  reference 
system.  The  estimated  observations  from  the  given  coordinates  x,  y,  z,  x,  y,  z 
are  obtained  as  follows: 

The  satellite  coordinates  in  a  topocentric  system  where  the  x-axis  points  east, 
y-axis  north  and  the  z-axis  completes  a  right-handed  system,  are: 


‘  xt' 

"  X  -  Xs  “ 

yT 

=  [s] 

y  -  ys 

ZT 

z  -  zs 

L  J 

and 

*T 

x+Wfc.y 

■  [B] 

y  -  WE  * 

• 

.ZT- 

• 

z 

(16) 


(17) 


where  xs,  ys,  zg  are  the  station  coordinates,  and  wE  is  the  rate  of  rotation  of  the 
Earth.  The  matrix  S  is  either 


or 


S  =  EM 


S  =  M 


depending  on  whether  the  transformation  is  performed  to  a  geodetic  or  a  geocen¬ 
tric  system. 

The  transformation  matrix,  M,  is: 


(U) 


-  sin(XGR 

-  \) 

cos  <AC,« 

-  A) 

0 

M  = 

-  C0S  <\,R 

-  \)  sin  4> 

-  sin  (Agr 

-  A)  sin  c|> 

cos  4> 

c°3 (Agr 

-  A)  cos  4> 

sin  (\CR 

-  A)  cos  tj> 

sin  4> 
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w 1 10 tv  A(  K,  \  and  <>  are  as  defined  previously. 

The  transformation  from  the  geocentric  to  the  geodetic  system  is  pci  lormcd 
by  the  matrix  F 

]  0  0 

E  j  o  l  -A4>  (in) 

i 

j^O  A<t>  1 

where  AO  *  -  6  and  is  a  small  angle. 

The  obtained  coordinates  x.r>  y ,,  zT,  xT,  yT  and  z.?  of  the  satellite  in  the  topocen- 

tric  system  allow  us  to  compute  the  observations  corresponding  to  the  estimated 
orbit. 

Observations  of  right  ascension  and  declination  are  usually  referred  to  a 
specific  celestial  system  (see  also  Section  VII  B)  defined  by  the  equinox  and 
equator  at  the  beginning  of  a  Bcsselian  year.  This  will  lie  discussed  in  more 
detail  in  Chapters  IV  and  VII.  At  the  present,  we  are  concerned  about  the 
transformation  of  the  satellite  position  from  the  basic  system  to  a  particular 
celestial  system.  This  can  be  done  by  utilizing  the  previously  given  precession 
and  nutation  matrices.  The  coordinates  of  the  satellite,  as  well  as  the  observing 
station  in  the  particular  celestial  system,  will  be: 


where  N  is  the  transformation  matrix  due  to  nutation,  P  is  the  transformation 

matrix  from  the  mean  equinox  of  date  to  that  of  1950  January  1 ,  ()h  UT  and  P 

is  the  transformation  matrix  to  the  mean  equinox  and  equator  of  the  celestial 

reference  system.  The  elapsed  days  from  1950  January  1,  0hUT  to  a  standard 
epoch  of  a  specific  Bcsselian  year  can  be  obtained  from  the  following  formula: 

d  -0,077  +  365.2422  (BY  -  1950.0)  (21) 

where  BY  is  the  Besselian  year;  1950.0  is  a  standard  designation  for  the  beginning 
of  the  Besselian  year  1950. 

If  the  celestial  reference  system  is  that  of  1950.0,  the  matrix  1J  need  not 

be  included  in  most  cases,  since  the  epoch  1950  January  1,  0hUT  is  very  close 
to  the  standard  epoch  1950.  0. 
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Another  method  is  to  obtain  the  right  ascension  and  declination  of  tne  satel¬ 
lite  in  the  basic  system  and  then  to  make  the  proper  corrections  for  nutation 
and  precession  to  obtain  the  right  ascension  and  declination  in  the  particular 
reference  system  of  a  standard  epoch. 
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IV.  SATELLITE  OBSERVATIONS 


A.  CLASSIFICATION  OF  OBSERVATIONS 

Satellite  observations  ean  be  obtained  generally  by  electronic  or  optical 
means.  The  electronic  systems  include  pulse  and  Doppler  radar  and  inter¬ 
ferometer-type  systems.  Optical  systems  include  visual  and  photographic 
observations.  Recently,  promising  results  have  been  obtained  by  laser  sys¬ 
tems  which  provide  angular  and  range  information. 

The  present  discussion  will  be  concerned  with  some  pertinent  aspects  of 
electronic  and  optical  systems  providing  observations  in  the  form  of  elevation, 
azimuth,  range,  range  rate,  and  right  ascension  and  declination.  In  addition, 
observations  of  elevation  and  azimuth  rates,  and  range  acce’e ration  will  be 
considered. 


B.  ELECTRONIC  OBSERVATIONS 

From  the  analyst's  point  of  view,  the  most  important  aspect  of  satellite 
observation  is  aeeuraey.  Errors  in  electronic  observations  arise  from 
several  sourees:  errors  in  the  electronic  and  meehanical  system;  errors  due 
to  refraction  and  aberration;  and  errors  in  the  physical  constants  of  the  station 
location.  Errors  in  the  electronic  system  are  due  to  causes  sueh  as  frequency 
drifts,  time  delays,  insufficient  resolution,  etc.  Other  errors  are  due  to  elec¬ 
tromechanical  systems.  These  errors  can  be  considerable.  Some  error  sourees 
eould  be  eliminated  to  a  large  extent  if,  instead  of  azimuth  and  elevation,  their 
rates  were  measured.  Thus,  the  precise  knowledge  of  the  true  meridian  and  the 
geodetie  vertical  is  not  critical.  The  bias  errors  in  the  rate  measurements 
should  be  smaller  than  in  angle  measurements. 

Range,  range  rate,  and  range  acceleration  measurements  rely  entirely 
on  the  electronic  system.  Thus,  a  source  of  some  highly  unpredictable  errors 
is  eliminated.  Moreover,  the  knowledge  of  the  true  north  and  vertical  is  not 
required.  Consequently,  these  measurements  can  be  mueh  more  accurate  than 
the  angle  measurements. 

One  of  the  ehief  sourees  of  error  in  electronic  measurements  is  refraction. 

It  is  generally  distinguished  between  tropospheric  and  ionospheric  refraetion. 

It  has  been  shown  (Ref.  G)  that  using  frequencies  in  the  kilomegacyele  range, 
the  ionospheric  refraction  is  practically  negligible  and  the  troposphere  refrae¬ 
tion  ean  be  ealeulated  with  fairly  good  aeeuraey,  particularly  if  the  loral  at¬ 
mospheric  conditions  are  taken  into  aceount. 

The  refraction  error  can  be  reduced  to  negligible  proportions  if  a  principle 
used  in  the  Transit  System  (Ref.  7)  is  utilized.  By  simultaneously  employing 
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two  coherent  frequencies,  the  actual  refraction  can  be  computed  to  a  lirst  order 
accuracy  and  the  measurement  corrected.  In  the  present  state  of  art,  the 
theoretical  accuracies  obtainable  for  range  and  range  rate  measurements  should 
be  of  the  order  of  5  m  and  0.2  m/sec.  However,  the  actual  accuracies  depend 
heavily  on  stringent  calibration,  maintenance  and  operation  procedures,  and,  of 
course,  the  type  of  equipment. 

C.  OPTICAL  OBSERVATIONS 

Some  of  the  most  accurate  satellite  observations  are  obtained  by  optical 
techniques.  Specifically ,  a  network  established  and  operated  by  the  Smith¬ 
sonian  Astrophysical  Observatory  (Ref.  8)  can  provide  observations  with  an 
accuracy  of  about  ±2  seconds  of  arc. 

Essentially,  the  precision  method  consists  of  photographing  satellites 
against  a  star  background.  Since  the  star  positions  can  be  determined  with  high 
accuracy  from  star  catalogs,  the  satellite  right  ascension  and  declination  can 
be  obtained  with  a  comparable  accuracy.  The  final  accuracy  of  optical  meas¬ 
urements,  however,  is  dependent  on  the  reduction  method  and  may  vary  by 
orders  of  magnitude.  Among  the  disadvantages  of  this  method  is  that  the  pre¬ 
cise  reduction  requires  elaborate  procedures  and,  therefore,  the  precise  data 
are  not  immediately  available.  Another  disadvantage  is  that  the  satellite 
can  be  photographed  only  at  certain  times  when  it  is  in  sunlight  and  the  observer 
in  darkness.  This  can  be  minimized,  if  the  satellite  carries  its  own  light  source 
or  by  laser  techniques.  One  of  the  chief  advantages  of  this  method  is  the  elimi¬ 
nation  of  many  sources  of  error.  The  chief  on-site  requirement  is  the  precise 
timing  of  the  exposure.  The  actual  computation  of  the  measurement  ean  be 
done  under  more  exacting  conditions. 


D.  OBSERVATION  CORRECTIONS 


The  radar  as  well  as  optical  observations  must  be  corrected  for  aberration 
and  refraction.  Because  of  the  relatively  high  velocities  of  satellites,  the  aber¬ 
ration  effect  will  be  significant  in  high  precision  measurements.  This  correc¬ 
tion  may  amount  to  a  few  seconds  of  arc  in  angular  measurements  and  several 
meters  in  range  measurements.  The  correction  is  done  on  the  basis  of  the 
satellite  velocity  vector  and  the  speed  of  light. 


A  much  more  significant  correction  is  required  for  retraction  effects.  In 
high  precision  measurements  of  right  ascension  and  declination,  obtained  by 
the  photographic  method,  the  major  part  of  the  refraction  is  corrected  indirectly. 
The  remaining  uncorrected  part  is  the  parallactic  refraction  which  arises  from 
the  fact  that  the  satellite  is  at  a  finite  distance,  while  the  stars  are  practically 
at  infinity.  Reference  8  gives  the  following  expression  for  parallactic  refrac¬ 
tion  correction : 


Afl  --435’.'0 


tan  z 


rT  cos  z 


-0. 1385  r 

e 


(22) 


20 


whore  /.  is  the  zenith  distance  aria  rT  the  range. 

The  refraction  correction  for  elevation,  range  and  range  rate  in  radar 
measurements  can  be  computed  by  the  following  method  (Refs.  9,  10,  11). 


The  troposphere  is  divided  into  m  incremental  layers  and  the  change  of  the 
index  of  refraction  in  3  layer  is  assumed  linear. 


The  average  index  of  refraction  as  proposed  by  the  NBS  Central  Radio 
Propagation  Laboratory  is  given  by  the  following  expression 

r?n  =  1  +  0.000313  exp  (-0.14386  1^)  (23) 

W'here  hnis  the  altitude  in  km. 

Designating 


N 


n 


-  l 


the  total  bending  through  the  troposphere  divided  into  m  incremental  layers  is 


y  = 


n=m 


I 

n=0 


2(Nn  -  Nn  +  1) 
tan  En4  tan  En^  <rad) 


(24) 
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The  elevation  angle  error  is  then 


AE 


ytan  E,„  -  (N0  -  Nm)  .  y2/2 
y  +  tan  Em  -  tan  E0 


(rad) 


and  the  range  error  is 

(N  +  N  )  (h  -  h  ) 

\  '  u  n  ♦  r  '  n  •  i  lv 

Z  Sin  En  +  sin  En  +  1 
n=0 


(l>5) 


(26) 


The  Doppler  velocity  can  be  corrected  to  a  first  order  approximation  by 

Ar  =  VAaTsin4>  (27) 

where  <i>  is  the  angle  between  the  line  of  sight  and  the  velocity  vector,  V. 

Aor  is  the  angle  between  the  line  of  sight  and  the  velocity  component  in  the 

direction  of  the  ray  path  at  the  target.  It  can  be  determined  from  the  expres¬ 
sion 


ik*  -  arc  cos 


cos  (EQ  -  AE) 


ai'c  cos 


C0SE- 


(28) 


The  above  expressions  are  valid  for  cases  where  the  frequency  is  in  the 
kiiomegacycle  range  and  the  ionospheric  refraction  effects  can  be  neglected. 
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V.  APPLICATION  OF  THE  MINIMUM 
VARIANCE  METHOD 


A.  FORMULATION  OF  THE  PROBLEM 

The1  knowledge  of  an  orbit  implies  the  knowledge  of  the  mathematical  model 
of  the  dynamical  system  and  certain  constants  associated  with  this  system. 

The  mathematical  model  can  be  developed  by  theory,  but  the  constants  must 
be  determined  by  experiment.  They  are  not  observed  directly,  but  can  be  de¬ 
termined  knowing  the  mathematical  relationships  between  the  observed  param¬ 
eters  and  the  constants.  In  our  analysis,  a  part  of  the  constants  is  represented 
by  the  so-called  state  variables,  whieh  in  our  ease  arc  the  instantaneous  orbital 
elements.  The  other  part  consists  of  the  various  physical  constants  of  the 
dynamical  and  the  observation  system.  Together  they  represent  a  generalized 
state  veetor. 

The  instantaneous  orbital  elements  themselves  are  not  constants.  However, 
they  are  rigorously  related  to  another  set  of  orbital  elements  at  some  other 
time  ehosen  at  an  epoch.  The  relationship  between  the  two  sets  is  defined  by 
the  mathematical  model.  Thus,  the  epochal  set  of  the  orbital  elements,  which 
are  constants  for  the  given  system,  completely  determine  the  instantaneous 
orbital  elements  at  any  other  time. 

An  approximate  orbit  can  be  obtained  by  assuming  a  simple  inverse- square 
central  force  field.  Various  methods  can  be  used  for  this  purpose  depending 
on  the  type  of  observations.  With  this  starting  orbit,  an  improved  estimate  of 
the  orbital  elements  and/or  the  various  physieal  constants  can  be  obtained  by 
more  sophisticated  methods.  The  application  of  the  Minimum  Variance  Method 
for  this  purpose  is  the  subjeet  ol  this  chapter. 

There  are  three  sources  of  errors  associated  with  the  orbit  estimation 
process:  incomplete  representation  of  the  dynamieal  system,  errors  in  the 
generalized  state  vector  and  errors  in  the  observations.  The  improvement  of 
the  dynamical  model  by  statistical  methods  is  beyond  the  scope  of  this  report. 
Statistical  knowledge  about  the  errors  in  the  generalized  state  vector  and  the 
observations,  however,  ean  be  utilized  in  orbit  improvement  from  observed 
data.  The  process,  which  is  called  statistical  filtering,  is  applied  to  obtain  a 
best  estimate  of  the  generalized  state  vector  on  the  basis  of  the  deviations  of 
the  actual  observations  from  the  estimated  orbit.  Once  a  good  estimate  of  the 
orbit  is  obtained,  past  or  future  instantaneous  orbital  elements  can  be  found  by 
prediction  and  smoothing  methods. 

Orbit  determination  has  always  been  a  major  problem  in  dynamieal  as¬ 
tronomy.  Methods,  notably,  the  Least  Squares,  have  been  in  use  for  approxi¬ 
mately  a  century.  With  the  advent  of  the  artificial  Earth  satellite,  these  methods 
were  adapted  for  the  new  application.  However,  problems  encountered  in 
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practical  space  engineering  applications  often  present  quite  different  situa¬ 
tions,  A  major  factor  is  the  electronic  computer. 

In  dynamical  astronomy,  relatively  few  observations  are  obtained  over  an 
extended  period  of  time.  In  many  engineering  applications,  a  large  amount  of 
observations  are  obtained  within  a  comparatively  short  time  period.  The  elec¬ 
tronic  computer  offers  immense  possibilities,  notably,  its  computing  speed 
and  storage  capacity,  However,  new  approaches  may  be  necessary  to  take  full 
advantage  of  these  capabilities  and  also  to  comply  with  some  of  the  numerical 
problems  arising  from  the  large  amount  of  data,  It  is  believed  that  the  Mini¬ 
mum  Variance  Method  offers  great  possibilities  in  many  space  engineering  ap¬ 
plications. 


B.  SOME  PRACTICAL  ASPECTS 

The  basic  difference  between  the  Minimum  Variance  Method  and  the  Least 
Squares  Method  is  that,  in  the  Minimum  Variance  Method,  the  orbit  is  continu¬ 
ously  updated  on  the  basis  of  each  new  observation  or  set  of  observations.  In 
the  Least  Squares  Method,  a  single  sotution  is  obtained  for  all  the  observations 
simultaneously.  In  each  case,  an  inversion  of  a  certain  order  matrix  is  in¬ 
volved.  The  order  of  the  matrix,  in  the  least  squares  case,  is  determined  by 
the  number  of  the  parameters  being  estimated.  On  the  other  hand,  the  order 
of  the  matrix  to  be  inverted  in  the  minimum  variance  case  corresponds  to  the 
lumber  of  simultaneous  observations.  ,cs  a  matter  of  fact,  all  the  simultaneous 
observations  need  not  be  processed  simultaneously,  and,  if  desired,  each  can 
be  handled  separately.  Thus  the  inversion  becomes  trivial.  Since  the  inversion 
of  high  order  matrices  is  not  a  simple  matter,  the  advantages  of  the  Minimum 
Variance  Method  are  obvious. 

Since  the  orbit  is  updated  on  each  new  set  of  observations,  in  many  applica¬ 
tions,  the  observations  need  not  be  stored  in  the  computer,  thus  saving  consider¬ 
able  storage  space.  The  updating  process  is  systematic  for  any  number  of  ob¬ 
servations,  either  too  few  or  too  many  to  obtain  a  deterministic  solution  for 
the  generalized  state  vector.  Also,  the  process  can  be  interrupted  at  any  time 
for  any  reason  giving  the  optimum  estimate  at  this  point,  and  resumed  at  a 
later  time.  The  Minimum  Variance  Method  allows  considerable  flexibility  in 
the  method  of  application  to  suit  the  particular  requirements. 

In  addition  to  the  important  practical  aspects,  the  method  is  able  to  consider 
correlated  measurement  errors  from  one  observation  time  to  another. 


C.  LINEAR  CONCEPTS 

The  theoretical  development  of  the  linear  filtering  theory  used  in  this  appli¬ 
cation  is  given  by  Kalman  in  Ref.  12  (see  also  Ref.  13).  An  application  of  the 
method  is  presented  in  Refs,  14  and  15,  The  following  presentation  will  deal 
with  the  special  case  of  satellite  orbit  determination  and  analysis. 


The  mathematical  model  of  the  dynamical  system  is  represented  by  second- 
order  differential  equations  of  the  form 


x(t)  -  f  [x(t),  x(t)] 


where 


are  the  three  position  coordinates. 

The  individual  terms  of  the  right-hand  side  will  be  considered  in  a  subse¬ 
quent  chapter.  For  the  time  being,  it  is  sufficient  to  know  that  the  equations 
are  nonlinear.  To  enable  the  application  of  the  linear  filtering  theory,  the 
equations  must  be  linearized  and  a  solution  obtained.  The  linearization  can  be 
accomplished  by  expanding  the  equations  in  a  Taylor  series  about  a  reference 
trajectory  and  retaining  only  the  first  order  terms.  Thus 


x(t)  =  f[xR(t),  xR(t)]  +  F(t)Ax(t) 


where 


afj  afj  afj 

w,  »  57  <•>  •  •  -H7  « 


F(t)=  3^(t) 


_^(t) 


.  ...  Ax„ 

Ax(t)=  3 


Ax  1 

r  3j 

The  differential  equations  for  the  reference  trajectory  are 
xR(t)  =  f  [xR(t),  xR(t)] 
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If  \vu  subtract  these  from  the  perturbed  equations,  we  obtain  the  linearized 
or  perturbation  equations 

Ax( t)  =  F(t>  Ax(t)  (32) 

The  three  second-order  equations  can  be  redueed  to  six  first-order  equations 
by  writing  them  in  a  standard  form. 

The  fundamental  solution  to  the  obtained  homogeneous  equations  is  called 
the  state  transition  matrix  and  designated  K  The  solution  can  be  written  in  the 
form 


Ax(t)  -  t(t,  to)Ax(t0)  (33) 

We  now  have  obtained  an  algebraic  matrix  equation  where  the  state  transition 
matrix  ?  relates  linearly  the  vector  Ax  at  time  t0  to  the  vector  Ax  at  time  t. 

in  a  general  case,  a  vector  Ax*  will  consist  of  all  the  parameters  to  be  deter¬ 
mined,  and  the  matrix  t*  will  be  of  a  corresponding  order. 

The  f*  matrix  possesses  some  important  properties,  which  are  summa¬ 
rized  as  follows: 


t*(t.  t)  =  I  =  unit  matrix 

(34) 

f*(tyU  *+(Vt,)  =**(tvt,) 

(35) 

?*'!(t2.t,)-  ?*(trt2) 

(36) 

in  practice,  the  matrix  t  for  the  six  orbital  elements  can  be  obtained  by 
several  methods.  In  our  case,  it  will  be  obtained  on  the  assumption  of  an 
unperturbed  orbit  in  an  inverse- square  central  force  field.  Experience  has 
shown  that  this  is  a  good  approximation  for  orbits  where  the  atmospheric 
and  gravitational  perturbations  are  not  significant.  For  very  low  altitude 
orbits,  the  approximation  may  be  insufficient  and  corrections  for  atmospheric 
effects  are  necessary  or  the  matrix  must  be  obtained  by  other  methods. 

The  equations  that  relate  the  observations  to  the  instantaneous  orbital 
elements  and  constants,  in  general,  are  also  nonlinear.  They  can  be  linearized 
by  a  similar  Taylor  series  expansion  about  a  reference  trajectory,  and  the 
equations  can  be  written  in  the  standard  matrix  form 

Ay(t)  =  H*(t)  Ax*(t)  (37) 

where  Ay  is  the  deviation  of  the  actual  observations  from  the  observations 
associated  with  the  reference  trajectory.  H*  is  the  matrix  of  partial  deriva¬ 
tives  of  the  observations  w'ith  respect  to  the  instantaneous  orbital  elements  and 
constants.  The  partial  derivatives  comprising  the  H*  matrix  can  be  obtained 
and  evaluated  at  the  observation  times  in  a  straightforward  manner. 
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The  third  important  matrix  used  in  the  filtering  equations  is  the  eovariance 
matrix.  At  any  point  in  the  filtering  proeess,  we  have  only  an  estimate  of  the 
instantaneous  orbital  elements  and  the  constants.  We  ean  eonsidcr  them  as 
random  sealars  or  components  of  a  random  veetor.  Thus  the  errors  in  the 
random  variables,  whieh  arc  deviations  from  an  expeeted  or  mean  value,  will 
also  be  random  variables  and  as  sueh  will  have  a  zero  mean. 


c  •  Ax*j  --  0 


The  estimate  of  the  individual  random  sealars  is  not  known  with  the  same 
aeeuracy.  However,  it  ean  be  described  by  the  so-ealled  variance  whieh  is 
defined  as 


o'  -  e  j[Ax*  -  eAx*r| 


where  e  Ax*  is  the  expected  value  of  the  deviations,  and  is  zero  in  our  ease. 
The  standard  deviation  is  defined  as  the  square  root  of  the  variance. 

The  individual  components  of  a  random  vector  can  be  affeeted  by  the  other 
components  whieh  is  ealled  correlation.  Thus,  instead  of  a  single  variance 
associated  with  a  random  sealar,  a  random  veetor  has,  in  general,  variances 
and  eovarianees.  The  eovarianee  matrix,  P*.  is  defined  as 

P*  -  cov|ax*,  Ax*  |  =  c|ax*  Ax*T|>  (40) 

T 

where  Ax*  means  the  transpose  of  Ax*. 

The  diagonal  elements  of  this  matrix  are  the  variances  and  the  off-diagonal 
elements  are  the  eovarianees.  If  the  components  of  the  random  vector  are 
uncorrelated,  the  off-diagonal  elements  will  be  zero.  If  the  eovarianee  matrix, 
P*  (to).  is  given  at  time,  tn,  the  eovariance  matrix,  P*  (t),  at  time,  t,  can  be 

obtained  by  use  of  the  state  transition  matrix  <p* 


P*(t)  =  eov|Ax*(t),  Ax*(t)|  =  e  |Ax*(t)  Ax*T(t)  j 


=  **(t,t0)P*(y  **‘(t,t0)  (41) 

Similar  considerations  apply  to  the  eovarianee  matrix  of  the  observation 
errors,  Q. 


D.  FILTERING  EQUATIONS 

The  statistical  filtering  theory  is  based  on  the  assumption  of  a  linear  multi¬ 
dimensional  dynamie  system,  whieh  ean  be  represented  by  the  following  model: 
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*'<t>  t(j)  x*(t,)  (42) 

y(t)  •  K*(t)  x*(t0)  ►  u(t)  (43) 

where  y*(i\  is  an  n-vector  and  represents  the  generalized  state  vector,  y(t)  is 
a  vector  representing  p  independent  measurements,  and  u(t)  is  a  p- vector 
representing  an  independent  gaussian  random  process  or  noise.  The  t>*  and 
H*  matrices  aiv  n  x  n,  and  p  x  n  matrices,  respectively,  which  have  been 
discussed  before.  Since  the  actual  system  is  nonlinear,  the  above  equations 
represent  the  linearized  or  perturbation  equations.  For  the  sake  of  simplicity, 
the  perturbations  are  represented  by  x  and  y  instead  of  Ax,  Ay. 


These  perturbations  are  referred  to  the  estimated  orbit.  Assuming  that 
at  time,  ,  an  estimate  ol  the  state  veetor  is  known  based  on  the  previous 

k  observations,  a  better  estimate  can  be  obtained  which  includes  the  k  +  1 
observation.  Under  the  linear  assumptions,  this  ean  be  expressed  as 

J  V  S*<V 

*  “'A,) [>'<W  -  »*<W  V  **<V]  («) 

whore  x  designates  the  estimate  of  x.  The  first  term  on  the  right-hand  side 
is  simply  the  transfer  of  the  state  vector  x*^)  at  time  *k  to  time  tkil  by  the  t* 

matrix.  Thus  it  represents  the  estimate  of  x*  at  tk  l  based  on  the  first  k 

observations.  The  second  term  represents  the  contribution  of  the  new  obser¬ 
vation  at  time  tk  ..  The  quantity  in  the  brackets  is  the  difference  between  the 

kth  observation  and  the  observation  based  on  the  estimated  orbit.  The  matrix 
K*(t.  )  is  a  weighting  matrix,  sometimes  ealled  the  optimum  gains  matrix, 

K*  1 

since  it  is  obtained  by  optimizing  a  loss  function. 


As  represented  in  the  above  equation,  the  assumed  model  of  the  dynamieal 
system  is  linear.  However,  in  orbital  analysis,  it  is  advantageous  to  use  the 
actual  dynamieal  system  represented  by  the  nonlinear  differential  equations 
of  motion  to  propagate  the  estimated  state  vector  from  one  observation  to  the 
next.  Also,  the  exaet  equations  can  be  used  to  obtain  the  observations  associ¬ 
ated  with  the  estimated  orbit.  With  these  modifications,  the  equation  reduces 
to 

X*<V.>  K*<V.)  [y<V,>  -  y'twj  <«> 

where  >'  is  the  actual  observation  and  y'  is  the  observation  based  on  the 
estimated  orbit.  Thus  the  linearized  equations  are  used  only  to  compi.te  the 
weighting  matrix  and  the  covariance  matrix.  The  advantages  of  this  approach 
are  that  the  process  is  less  sensitive  to  errors  in  the  initial  conditions.  The 
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nominal  orbit,  being  continuously  updated,  does  not  deviate  excessively  from 
the  true  orbit.  Consequently,  the  chance  of  violating  the  linearity-  assumptions 
is  minimized. 

Kalman  (Kef,  12)  obtains  an  optimal  weighting  matrix  K*  for  an  independent 
gaussian  random  process  utilizing  the  Schmidt  orthogonalization  procedure  in 
a  multidimensional  space.  The  optimal  filter  or  weighting  matrix  is  obtained 
utilizing  the  first  and  second-order  statistics:  the  expectation  of  the  state  vec¬ 
tor  and  the  covariance  matrix. 

-i 

K‘<U  *  P*<W  H’X.,)  [H*(tk.,>  P‘<‘k,>  »*'\„>]  <«> 

where  P*  is  the  covariance  matrix  of  the  state  vector  and  H*  is  a  matrix  of 
the  partial  derivatives  of  the  observations  with  respect  to  the  components  of 
the  state  vector.  The  matrix  P*(tkl)  is  obtained  transferring  the  P*’^)  matrix 

by  means  of  the  state  transition  matrix  ?*  and  adding  the  error  matrix  Q* 

p*<‘k.,>  -  V  p*\>  V  -  Q*  <tw>  (47) 

Knowing  the  weighting  matrix  K*(t  ),  a  new  covariance  matrix  at  time  t  , 

K-l  k*l 

which  includes  the  statistics  of  the  new  observation,  can  be  obtained  by  the 
recursion  equation 

-  P*(tk  j)  -  K*(tk  j)  H*(tk.j)  P*{tk+,)  (48) 

Since  the  state  vector,  generally,  can  be  composed  of  many  parameters,  the 
matrices  will  be  of  a  large  order,  affecting  the  numerical  operations. 

In  Ref.  14,  it  is  shown  that  a  simplification  of  the  matrix  equations  can  be 
achieve^  if  the  errors  in  the  observations  are  uncorrelated  from  one  observa¬ 
tion  time  to  another.  It  is  likely  that  for  the  same  station,  the  observations  will 
be  correlated  in  this  sense  because  of  imperfect  calibration,  etc.  In  most 
cases,  it  will  be  sufficient  to  represent  such  errors  by  an  algebraic  function. 
Then  the  constants  in  this  function  can  be  regarded  as  bias  errors,  included 
in  the  state  vector,  and  ihus  estimated  along  with  the  other  constants.  The 
remainder  of  the  errors  then  can  be  considered  as  a  Gaussian  random  vector 
and  thus  uncorrelated  between  two  observation  times.  With  this  assumption, 
the  equations  reduce  to 

*<‘M>  -  K<‘k„>  [y'lW-y’tt,.,)]  (49) 

K(‘k.,)  ‘  puk. ,)  HT(tk.,)[H(tk.,)  Pit,.,)  HT (tk„)  4.  Q(tk„)]' '  (50) 

<5l> 


2!) 


hi  )  P(t 


K-l' 


Kit,.,)  H( ihil)  Pit 


Lkq' 


(hd) 


where  all  the  matrices,  except  the  Q  mat  x,  pertain  to  the  system  state  vector 
only  .  1  he  Q  matrix  represents  the  covariance  matrix  of  the  observations. 

II  the  slate  vector  consist.-  of  the  six  orbital  elements,  the  matrices  will  never 
be  largei  man  ii  x  <>.  The  matrix  to  be  inverted  will  is.  of  an  order  equal  to 
the  numbir  of  the  simultaneous  observations. 

The  inversion  of  the  matrix  is  an  important  operation  and  deserves  further 
analysis.  It  the  observations  are  linearly  independent,  i,  e.  ,  if  none  oi  the 
observed  scalar  random  variables  is  a  linear  combination  of  the  others,  the 
matrix  will  be  invertible  whenever  P  and  Q  are  positive  definite.  By  virture 
of  definition,  the  covariance  matrices  P  and  Q  are  positive  definite,  and  thus 
the  combination  is  invertible.  The  practical  aspects  of  inverting  a  matrix, 
however,  are  quite  different  front  those  of  pure  theoretical  considerations. 

H  lias  1  ieeii  pointed  out  in  Kef.  l(i  that  the  differential  correction  matrix  tends 
to  become  singular  as  the  time  are  of  the  filtering  process  increases,  The 
singularity  is  apparently  also  affected  by  the  choice  of  the  orbital  elements. 

In  the  Least  Squares  Method,  this  will  prevent  a  solution.  In  the  Minimum 
Variance  Method,  the  situation  is  somewhat  different.  First,  as  pointed  out 
earlier,  tin  order  of  the  matrix  to  hr  imeried  is  c-ejuaj  in  die  number  of 
simultaneous  observations  and,  if  desired,  all  the  observations  at  a  particular 
lime  iv ri  not  be  considered  simultaneously.  The  inversion  thus  can  he  made 
trivia!.  Secondly,  in  mod  practical  eases,  the  observations  can  be,  indeed, 
i  wnsidi  ri  d  independent ;.  .  i  > v ■■  •n.c.lntec!  from  one  ohsorv  bon  time  to  another. 
The  v  ioainx  thus  will  lie  .1  diagonal  m:itr,x.  ic.pi  serai  ii,’,  tie  va-Munces  of  the 
ob'-ei  .  tit  <i.s..  As  the  filu.jo.ig  progresses,  the  orbit  wall  become  Known  with 
higher  accuracy.  fn  other  words,  the  eteinenhs  of  matrix  P  will  assume  smaller 
men  :v  ii  \  .dues.  Co».->cquemi\  ,  the  dominance  of  the  Q  manic;  w  ill  become 
more  i  .  r -minced,  and  since  it  1-  a  diagonal  matrix,  no  inversion  problem  will 
be  ein  a  note  red.  Thus  the  .amice  of  the  orbital  element.-,  is  not  t  1  itieal  and 
may  be  made  on  the  basis  oi  otner  considerations,  in  die  present  case,  the 
orbital  elements  are  the  position  and  velocity  coordinates.  In  a  program 
pm  scmlc  under  development .  no  inversion  diifieuities  have  been  encountered, 


K.  it!:..)  liCTtoN  OF  OB.SF.RVATIONS 

Tin  .  b  litmiin  Variane..'  Method  is  applicable  to  normally  distributed  random 
■  i  101’  .  it  is  inevitahh  liiat,  in  tin;  observation  data,  liiere  will  be  observations 
v.oirl.  t  <r  n:  reason  or  e.otu.r  will  he  greatly  in  error.  As  seen  they  should 
ml  be  hided  in  the  estimation,  anti  some  kind  of  criterion  should  be  employed. 
Ills  ;  be  lone  as  follows. 

At  art  lime  during  the  filtering,  an  observation  with  a  standard  error  must 
•  ail  wilim;  a  region  represented  by  the  sum  of  the  covariance  matrices  of  the 
estimated  orbit  and  the  observations. 


;to 


R 


f»3) 


H(tk)  P(tkj  H1  (tk)  +  Q(tk) 

This  covariance  matrix  represents  an  error  ellipsoid  of  the  observations 
(Ref.  17).  Because  of  the  covariance  elements,  the  principal  axes  of  this 
ellipsoid  do  not  coincide  with  the  axes  system  in  which  the  observations  are 
represented.  Since  the  matrix  is  real  and  symmetric,  it  can  be  diagonalized 
by  a  similarity  transformation.  The  diagonal  elements  thus  would  represent 
the  excursions  of  the  observations  in  the  direction  of  the  principal  axes.  How¬ 
ever,  the  identity  of  the  original  observations  would  be  lost. 

Leaving  the  R  matrix  intact,  we  can  partition  it  to  a  single  element,  and 
consider  only  the  diagonal  elements.  These  arc  called  marginal  deviations, 
and  represent  a  case  where  the  deviation  of  an  individual  observation  is  con¬ 
sidered  under  the  assumption  that  all  the  other  deviations  can  be  infinity.  In  a 
general  case,  instead  of  a  hyperellipsoid,  this  will  give  a  rectangular  hyper- 
parallclopiped.  The  square  roots  of  the  diagonal  elements  of  matrix  R  thus 
will  represent  the  standard  deviations,  q,  under  these  assumptions.  Now,  any 
actual  observation  which  exceeds  na,  referred  to  the  estimated  orbit,  can  be 
rejected.  In  the  case  of  a  good  representation  of  the  mathematical  model  and 
knowledge  of  the  observation  errors,  n  could  be  3  (3o).  It  must  be  emphasized 
that  the  covariance  matrices  do  not  represent  absolute  numbers  and,  therefore, 
should  be  tieated  accordingly. 


F.  DETERMINATION  OF  CONSTANTS,  TYPE  1 

The  previously  given  matrix  equations  can  involve  operations  with  large 
order  matrices,  in  case  constants  are  estimated  simultaneously  with  the  six 
orbital  elements.  However,  recognizing  the  nature  of  the  constants,  certain 
simplifications  can  be  introduced  and  the  order  of  the  matrices  reduced. 


There  is  a  certain  class  of  constants  which  are  not  functions  of  certain 
observations.  In  this  case,  the  partial  derivatives  of  the  observations  with 
respect  to  the  constants  are  all  zero.  Constants  of  this  type  include  the  coef¬ 
ficients  of  the  zonal  harmonics,  exospheric  temperature  of  the  atmosphere, 
Earth's  mass,  etc. ,  in  combination  with  observations  of  right  ascension, 
declination,  azimuth,  elevation,  range,  range  rate,  azimuth  rate  and  elevation 
rate. 


Returning  to  the  previously  derived  equations  for  x+,  K*,  and  P*,  we  see 
that  no  difficulty  is  encountered  writing  ic*  in  a  partitioned  form 


(54) 
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when.1  e,  n,  now  denote  estimates  oi  perturbation  vectors  representing  the 
six  orbital  elements,  constants  and  measurements,  respectively.  Thus  the 
veetors  can  he  eslimaud  mil.  ;  set  a  tenth .  It  remains  to  investigate  the  K*  and 
P*  matrices.  We  can  write 
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I 

where  P*(t  )  is  in  a  partitioned  form 
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P*(tk) 


r  p.  p.  p. 
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P1  P:  P' 

\n  cn  n 


{55) 


(56) 


and  Px,  P1,  Pn’ai'e  the  covariance  matrices  of  tin1  orbital  elements,  constants, 

and  observations,  respectively,  The  submatrices  with  double  subscripts 
designate  the  covariances  between  the  throe  groups.  Since  a  covariance  matrix 
must  be  symmetric,  the  nth  dnyonal  subma trier  s  must  be  transposes  of  each 
other.  The  state  transition  matrix  f*(t  I  )  can  be  expressed  as 

t\  •  I  r  Y 


i  f* 

?  0 

‘x 

xc  1 

:*(t. 

..  t) 

1  0 

I  0  1 
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■  0 
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(57) 


where  the  ?\  matrix  is  tin*  state  transition  submatrix  of  the  orbital  elements. 
The  state  transition  sulnnat n.  o!  the  constants  is  a  unit  matrix  I.  fcx,is  a  sub- 

matrix  relating  the  state  of  the  elements  at  time  t  ,  as  affected  by  the  state 

of  the  constants  at  t^.  The  other  submutricos  are  zero  because:  (1)  there  is 

no  correlation  between  the  observation  errors  at  time  tk,  which  are  assumed 

Gaussian  random  error1-,  and  the  orbital  elements  or  constants  at  time  tk 

(2)  the  constants,  arc  not  nlfeelcd  by  errors  in  the  orbital  elements’,  and  (3)  the 
observation  errors  are  uncorrelated  from  one  observation  time  to  another. 
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and  represents  the  observation  errors. 


Now,  performing  the  required  operations  we  obtain 

[<*x  K+  «x<0  *J  +  <»xP*e  +  *xcpc  >  flc]  [»xPU  •  V? 
[pT  tT+  p'  tT  1 

[xc  X  c  xcj 


p*(tk.^ 


xc 
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P  P  0 

X  xc 

P  7  P  0 

xc  c 

0  0  Q’ 
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Thus  the  transfer  of  the  covariance  matrix  P*'(t  )  can  be  accomplished  in 


parts 


Pxc  "  KKc  *  *xcPc 
Px  =  <«'  +  fx,P»T)  +  P«c  * 


p  =  P' 

c  c 


(60) 

(61) 

(62) 


Note  that  Pc  =  Pc. 

In  the  expression  for  K*(tk+1),  the  quantity  in  brackets  is 
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H*  P*  H*T  =  H„  P„  HP  +  H„Q'H^ 


In  our  ease  Hc  =  0,  and  performing  the  matrix  multiplication 

R 

1 1  '  ii 

Thus  the  weighting  matrix  is 
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or 


K*  =  P*  H*T  R*1  = 
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(63) 

(64) 
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Performing  the  operations  in  the  recursion  equation  we  obtain 

p*>  =  p*  _  K*  P*  = 
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Now  the  updating  of  the  P*  matrix  can  be  done  by  parts.  Thus  the  covariance 
matrix  of  the  orbital  elements  is 


P’ 


=  p  -  p  h  1  r  1  h  p  =  p  -k  h  p 

X  XX  XX  X  XXX 


(65) 


The  updating  of  the  covariance  matrix  of  the  constants  is  accomplished  by 
the  equation 


P*  -  P  -  P  T  H 1  R* 1  H  P  =  P  -  K  H  P 

c  r  xr  x  x  xc  *-  c  x  xc 


(66) 


It  should  be  noted  that  the  second  term  on  the  right-hand  side  is  strongly 
dependent  on  the  correlation  between  the  constants  and  the  orbital  elements. 

If  the  correlation  is  weak,  the  updating  will  be  very  ineffective,  as  will  be  the 
estimation  of  the  constants. 


Similarly,  the  covariance  submatrix  P  will  be  updated  as  follows 


P'  ■•=  P  -P  H  R‘HP  =  P  -KHP 

XL  XC  X  X  X  XC  XC  XX  XL 


(67) 


The  other  submatrices  P  and  P  need  not  be  updated  because  they  will 

xn  cn  r  J 

vanish  in  the  transformation. 


G.  DETERMINATION  OF  CONSTANTS,  TYPE  2 

A  similar  simplification  of  the  filtering  equations  can  be  accomplished  if 
the  constants  to  be  estimated  are  not  related  to  the  orbital  elements,  i.  e.  , 
they  are  not  included  in  the  equations  of  motion.  Constants  of  this  type  include 
the  coordinates  of  observing  stations  and  bias  errors  in  the  observations. 
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Annin  wc  can  write  the  covariance  matrix  in  a  partitioned  form: 


P*’(tk)  = 
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where  the  submatrices  are  defined  previously. 


The  state  transition  matrix  now  will  be 


®x  °  0 
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I  0 

0  °_ 


(69) 


Regarding  the  correlation  from  one  observation  time  to  another  we  have 
made  the  same  assumptions  as  in  the  previous  case  with  the  following  modifi¬ 
cation.  The  observation  errors  now  are  assumed  correlated  and  represented 
by  an  error  function.  The  constants  in  this  function  are  represented  in  the 
covariance  matrix  Pr  and  estimated  together  with  the  other  constants  of  this 

type.  The  uncompensated  observation  errors  arc  considered  Gaussian  random 
errors  and  as  such  uncorrelated  from  one  observation  time  to  another. 

The  transition  submatrix  $  will  be  zero  because  the  constants  do  not  enter 

into  the  equations  of  motion.  The  covariance  matrix  Q*  is  the  same  as  pre¬ 
viously. 


Performing  the  required  operations  we  obtain 


(70) 


and  the  transfer  of  the  covariance  matrix  can  be  done  in  parts 


u 

(71) 

=  ?  P’ 
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x  xc 

=  Pc 

(73) 

The  weighting  matrix  K*  can  be  obtained,  first,  considering  the  expression  in 
the  brackets. 
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II*  P*  H*r  -  (H,  P,  •  Hcpj  )Hr+  (Hx  pxc  *  Hcp(.  *  H„Q'  H„r 


R 


The  weighting  matrix  is  then 
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K>  =  ft  *PschJ  ;r-‘ 
Kt=(Pj  Hi  +  PC  Hl)R-: 


(74) 

(75) 


The  recursion  equation  then  can  be  obtained  as  follows: 
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Which  gives  after  performing  the  required  operations 

K  =  p*  -  Kx(HxPx  *HCPXJ)  (76) 

P’c  -  Pxc-  KX(HX  Pxc  +  HcPr  )  (77) 

P;  =Pr  Kc(Hx  Pxc.  +  HCPC  )  {78) 


Thus  the  order  of  the  matrices  is  reduced,  and,  for  the.  diagonal  submatrices, 
it  will  be  equal  to  the  number  of  the  orbital  elements  or  constants,  respectively. 
For  the  covariance  matrix  P*,  the  submatrices  forming  the  rows  and  columns 
are  simply  transposes  of  each  other.  The  submatrices  I*'  and  I*’  and  their 

transposes  in  P*'  need  not  be  computed  since  they  will  vanish  in  the  transfor¬ 
mation.  In  either  case,  the  maximum  order  of  matrix  R,  which  must  be  in¬ 
verted,  is  equal  to  the  number  of  simultaneous  observations. 

The  equations  as  derived  in  this  and  the  previous  section  are  for  a  simul¬ 
taneous  estimation  of  the  orbital  elements  and  constaits.  Obviously,  the  whole 
process  can  be  separated  in  two  parts.  First,  the  constants  may  be  assumed 
known  and  a  best  fitting  orbit  determined  as  in  a  normal  orbit  determination 
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routine.  The  residuals  then  can  be  attributed  to  the  constants  and  the  filtering 
repeated  considering  only  the  equations  pertaining  to  the  constants.  The  corre¬ 
lation  between  the  orbital  elements  and  the  constants  thus  will  be  ignored  but 
the  estimation  process  will  be  simplified. 
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VI.  STATE  TRANSITION  MATRIX 


A.  FORMULATION  OF  THE  PROBLEM 

One  of  the  matrices  in  the  filtering  equations  is  the  state  transition  matrix 
of  the  six  orbital  elements.  This  matrix  can  be  obtained  for  the  exact  mathe¬ 
matical  model  by  the  so-called  secant  technique,  i.e. ,  perturbing  ihc  elements 
one  at  a  time  and  obtaining  a  solution  by  integration.  For  observations  close 
together,  and  using  the  Minimum  Variance  Method  in  a  point-to-point  mode  as 
outlined  previously,  this  can  be  a  very  efficient  method.  If  the  observations 
are  far  apart,  the  method  is  very  time  consuming. 

The  matrices  used  in  the  filtering  equations  are  not  arrays  of  absolute  num¬ 
bers,  and  so  do  not  require  absolute  precision.  Therefore,  a  good  approxi¬ 
mation  of  the  actual  dynamical  system  is  permissible.  Experience  has  shown 
that  such  good  approximation  is  a  simple  inverse-square  central  force  field. 
The  resulting  orbits  are  a  circle,  ellipse,  parabola  or  hyperbola, depending 
on  the  eccentricity.  Even  with  this  approximation,  the  analytical  solution  has 
presented  considerable  challenge.  As  a  result,  many  analytical  methods  have 
been  developed  and  published.  A  pure,  closed  form  analytical  solution,  how¬ 
ever,  is  not  always  the  most  satisfactory  for  the  electronic  computer.  The 
computer  is  most  efficient  for  repeated  solutions  of  simple  arithmetic  equa¬ 
tions,  which  save  storage  space  and  computation  time.  The  method  that 
follows  has  been  developed  with  these  considerations  in  mind. 

The  problem  can  be  stated  as  follows.  Given  the  six  orbital  elements  in 
the  form  of  rectangular  coordinates  x^  yp  zv  x,,  yp  iv  at  time  tp  find 
a  state  transition  matrix  which  is  defined  as  one.  relating  small  perturbations 
of  the  state  at  time  t  (  to  the  resulting  perturbations  at  time  t.,.  A  fundamental 
solution  will  be  obtained,  first,  for  an  elliptic  orbit,  and  then” extended  to 
circular,  parabolic  and  hyperbolic  orbits.  Since  orbits  in  a  central  force 
field  are  planar,  the  solution  can  be  obtained  in  three  steps:  ^1)  in-plane 
perturbations:  (2)  out-of-plane  perturbations:  (3)  transformation  to  the 
original  axes  system. 


B.  ELLIPTIC  CRBITS 

The  direction  cosines  of  an  axes  system  in  which  the  xu-axis  is  directed 
toward  the  point  on  the  orbit  at  time  t  ,  yu-axis  is  in  the  orbital  plane  such 

thaty  >  0,  and  z^-axis  completes  a  right  hand  system  are,  first  for  the  x  - 
axis 


3‘J 


(79a,  b,  c) 


The  direction  cosines  for  the  z -axis  are  obtained  by  taking  the  vector  product 
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Similarly,  the  direction  cosines  of  the  y -axis  are  obtained  from  the 
vector  product  of  the  unit  vectors  in  the  direction  of  zu-and  x^-axes. 
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The  velocity  components  in  the  new  planar  axes  system  arc 
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also  y  =z  -  z  -  0 
i  *1  “l 


(82) 
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We  can  now  compute  a  set  of  orbital  elements  which  define  the  planar  orbit 
(see  Appendix) . 
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A  and  B  give  eccentricity 
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Rotating  the  planar  axes  system  so  that  the  x-axis  points  toward  the  peri- 
center,  the  direction  cosines  of  the  new  system  are 


£ 


3 


A 

e 


(88a,  b) 


and  the  normalized  coordinates  of  the  initial  point  in  this  axes  system  are 
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(89a,  b) 


We  can  now  find  the  complementary  angle 
e  +  (l  -  e2)  x 


9, =  arctan 


(i  -  W2  y, 


The  eccentric  anomaly  is 

Ei  =  T  ~  Qi  if  h  >  0 
E,  =  I  *  -  9,  if  y,  <  0 


(90) 


(91a) 

(91b) 


The  mean  anomaly  for  the  initial  point  is  from  Kepler's  equation 

M,  =  E.  -  e  sin  E,  (92) 
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Now  the  eeccntric  anomaly  for  the  final  point  can  be  obtained  by  an  iteration 
method.  For  a  first  approximation 
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(93) 


where 
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and  At  =  t.,  -  tx 


and  the  solution  ean  be  obtained  to  the  desired  degree  of  accuracy  by  success¬ 
ively  computing 
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At  -  Atp 

E,  = 
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where  the  subscript  o  designates  the  previous  estimate. 

The  position  and  velocity  of  the  final  point  in  the  planar  axes  system 
ean  be  obtained  by  solving  the  following  set  of  equations 
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This  process  must  be  repeated  either  four  or  eight  times,  perturbing 
successively  x^  ,  ,  x^,  ,  and  starting  with  Eq.  (84).  Since  the 

perturbation  of  means  the  reorientation  of  the  axes  system,  the  perturbed 
•  * 

velocities  x  and  v  must  be  transformed  in  the  new  axes  system,  before 
solving  for  the  elements  A,  B,  C,  by  the  matrix 
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where  Ay  is  the  perturbation  and 
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(101) 
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After  obtaining  the  final  values  x  ,  y^  ,  xw  ,  yw  for  this  ease,  they  must 

2  2  2  -J 

be  transformed  into  the  original  system  by  multiplying  by  the  transpose  of  this 
matrix. 


Each  perturbation  will  give  the  partials  of  x^ 

^2 

to  the  particular  perturbation.  Thus 
Ax  x  -  x 


y  ,  x  y  with  respect 

u2  2  ~2 


0 
1  1 


where  x 


2n 


Ax 


A  x 


etc, 


*i  l 

is  the  nominal  value  of  x^  . 

U2 


A  state  transition  matrix  for  a  planar  orbit  in  the  x_  -  y^  -  axes 
system  can  be  written  as 
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<t> 

4a 

0 

d1 

*52 

0 

*-a 

o4 

<> 

3  5 

0 

0 

L. 

0 

6  , 

0 

0 

*«6_ 

(102) 


The  first,  second,  fourth,  and  fifth  columns  of  the  matrix  have  been  obtained 
by  the  perturbation  technique.  The  indicated  elements  in  these  columns  n  it  be 
zero  because,  for  a  planar  orbit,  the  in-plane  perturbations  cannot  cause  out- 
of-plane  deviations.  It  remains  to  determine  the  third  and  the  sixth  column 
which  is  due  to  the  out-of-plane  or  Az  and  Az^  perturbations.  The  mam 

effect  of  these  perturhations  is  to  tilt  the  resulting  orbit  with  respect  to  the 
original  orhital  plane.  For  small  perturbations  of  z^  ,  and  zw  ,  the  tilting 

angle  and  the  increment  in  total  velocity  is  small.  To  a  first  order  approxi¬ 
mation,  assuming  that  the  cosine  of  small  angles  is  equal  to  one,  they  will 
have  no  effect  on  the  xw  and  y^  coordinates.  Thus  the  indicated  elements 
in  columns  three  and  six  are  assumed  zero. 
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Considering  a  new  axes  system  with  the  x  -axis  displaced  by  A/:  and 

the  new  orbital  plane  determined  by  the  new  velocity  vector ~v  *Az 

we  can  obtain  the  direction  cosines  of  the  new  z'_-axis  neglecting  higher 
order  terms 


i  i 


Az - — 

y. 


K  =  i 

Because  of  orthogonality,  a  projection  of  coordinates  in  the  new  system 
on  the  original  z  -axis  will  be 


x  + 
"2 


x  Az 

w  <~ 

_1 _ 1 

x  y 


y  -  z 


We  previously  established  that,  because  of  the  small  angles,  £  =  x  and 

^2  ~2 

y  =  y  to  a  first  or^^r  approximation.  Now  if  z  =0  (for  a  planar  orbit) 


AZ  = 


x.  y. 

_J _ 2_ 

x  y 
wi  l 


Az  +  — 

ui  y. 


By  comparing  this  equation  to  the  third  row  in  matrix  <t>,  we  obtain 


4>  - - 

33  Az 


x  y 

1 _ 2_ 

X  V 

“I  W1 


(103) 


6  = - 

3t>  Az 


(104) 
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By  :t  similar  process  we  obtain 


Az  x  x  y 

^fi3_Az  ~x  yw 

J  U1  wi  wi 

Az  y 

w2  2 

titi  —  1 — 


(105) 


(106) 


Thus  all  the  elements  in  matrix  4>  are  known.  Forming  a  matrix  X  from  the 
previously  computed  direction  cosines 


5,  ?2  ?3  0  0  0 

r).  r]  n  o  0  0 

C,  ?2  ?3  0  0  0 

o  o  o  s,  ?2  ?3 
o  o  o  vi,  n2  n3 

o  o  o  ?,  ?2  C3 


(107) 


We  can  now  obtain  the  state  transition  matrix  in  the  original  rectangular  coordi¬ 
nate  system 


t  =  X  ❖  Xr  (108) 

The  general  procedure  is  valid  for  all  eccentricities  with  exception  of  the 
particular  computations  as  outlined  next. 

C.  CIRCULAR  ORBITS 

In  the  expressions  for  the  direction  cosines  ?a  and  r)a,  the  eccentricity  e 

j  j2 

appears  in  the  denominator.  However,  because  e  is  obtained  from  (A2  +  B2)  , 

no  numerical  difficulties  will  be  encountered  for  small  eccentricities  as  long  as 
a  sufficient  number  of  significant  figures  are  carried.  Obviously,  circularity  is 
a  relative  matter.  Thus  at  some  point  the  orbit  can  be  assumed  circular. 

Since  a  circular  orbit  has  no  pericenter,  we  can  assume  an  orientation  of  the 
x-axis  to  coincide  with  the  initial  point.  Therefore,  the  direction  cosines  are 

5a=-i,  n.=  o 

and  the  final  eccentric  anomaly  is 


(109a,  b) 


(110) 
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The  rest  ol'  the  solution  is  the  same  as  for  the  elliptic  case, 

D.  PARABOLIC  ORBITS 

A  parabolic  orbit  has  an  eccentricity,  e  =  1.  The  chances  of  a  pure  parabolic 
or  a  pure  circular  orbit  occurring  in  computations  are  low  and,  as  pointed  out 
previously,  in  fact,  is  a  relative  matter.  In  the  assumed  rectangular  coordinate 
system,  the  solution  for  a  parabolic  orbit  can  be  obtained  in  a  closed  form. 

First,  the  normalized  area  swept  out  by  a  radius  vector  from  pericentcr  to  the 
initial  point  is  (see  Appendix) 

y 3  y  x 

J  i  i 

*.=i — —  <iu> 


and  the  area  swept  out  to  the  final  point  is 

A,  =  V  At_ 
Kc 

where 


from  which  the  coordinates  of  the  final  point  are 


/  o  \  (i 

r2T/3  r  /  2  \l/ 

2"jl/3 

6A,  h  (36A“  •»  11 

+  6A2  -  (36A;2  4  lj 

The  rest  of  the  solution  is  the  same  as  for  the  elliptic  case. 

E.  HYPERBOLIC  ORBITS 


(112) 


(113a) 

(113b) 


Hyperbolic  and  elliptic  orbits  are  the  most  important  orbits.  A  hyperbolic 
orbit  is  one  with  eccentricity,  e  >  1.  The  solution  for  a  hyperbolic  orbit  must 
be  obtained  by  iteration.  We  can  write  the  normalized  area  swept  out  by  the 
radius  vector  from  pericenter  to  the  initial  point 


Ai  = 


ey, 


(e -1) 


1/2 


In  B. 


(114) 


where 


B,  -  e2  -  (e2-x)  xi  -  (e2-l)U2  y{ 


For  a  first  estimate  of  the  final  point  we  assume  a  parabolic  orbit.  The  esti¬ 
mate  of  the  swept -out  area  to  the  final  point  is 


A,  -  Ai  - 


At 


whe  rc 


K  .  =  —  (--) 
^  e2-l 


2(e'-l)Kt. 

3v  1/2 


(115) 


The  estimated  y.,  and  x2  are  obtained  from  Eqs.  (113a) ,  (113b)  and  an  improved 
A,,  obtained  from 

a,  -  ey,  +  -2— a  In  B,  (116) 

'  (e2-l) 


whe  re 

b2  =  e  -  (e2-1)*,  -  (e2-l)l/2  y2 


The  time  corresponding  to  the  estimated  final  point  is 

(A2  ~Ai) 


At0  =  Kc 


(117) 


The  iteration  is  done  on  y^  successively  solving  Eqs.  (116).  (117),  (118),  and 
(119)  until  the  desired  degree  of  accuracy  is  reached. 


y, =  X,  4 
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2  2 

1  Me  -l)y2 


e2-l 


(118) 

(119) 


The  remaining  solution  is  the  same  as  for  the  elliptic  case. 

Experience  has  shown  that  the  initial  estimates  of  x2  and  y2  in  the  hyperbolic 

case  are  unimportant  and  even  grossly  inaccurate  estimates  will  give  rapid 
convergenee  using  the  above  method.  Therefore,  the  use  of  the  parabolic 
solution  for  the  initial  estimates  is  justified. 
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VII.  EQUATIONS  OF  OBSERVATIONS 


A.  REQUIREMENTS  FOR  ANALYSIS 

To  obtain  the  residuals  or  differences  between  the  observed  and  estimated 
observations,  it  is  necessary  to  compute  the  associated  observations  from  the 
estimated  orbit.  In  addition,  the  filtering  equations  require  the  H-matrix  which 
has  been  defined  as  a  matrix  of  the  partial  derivatives  of  tie  observations  with 
respect  to  the  state  variables.  As  shown  previously,  the  h-matrix  can  be  par¬ 
titioned  to  simplify  the  matrix  operations.  It  was  also  noted  that  for  a  large 
class  of  constants  the  Hr -submatrix  was  either  zero  or  a  unit  matrix.  The 

relationship  between  the  observations  and  the  orbital  elements  as  well  as  the 
associated  Hx -matrix  remains  to  be  determined. 


B.  ANALYTIC  EXPRESSIONS 


To  obtain  the  expressions  relating  elevation  and  azimuth  to  the  orbital  ele¬ 
ments,  we  will  utilize  the  equations  given  in  Section  III-F.  Knowing  the  satellite 
coordinates  xT,  y  ,  zT  in  the  topocentric  axes  system  as  defined  previously, 

elevation  and  azimuth  can  be  expressed  as 


E  = 


A 


=  arc  tan 


(120) 


(121) 


and  the  partial  derivatives  with  respect  to  the  six  orbital  elements  in  the  form 
of  position  and  velocity  coordinates  are 


9E  _  9E  9E  = 
dx  9y  9  i 


(122) 


(123) 
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i~ 


and 


(124) 


a  a  a  a  ...  a  a 

ax  3y  3z 

where 


(125) 


and  the  matrix  S  has  been  defined  previously  in  Section  III-F. 

The  expressions  for  the  estimated  right  ascension,  o,  and  declination,  <5 , 
must  be  obtained  in  the  particular  celestial  system  in  which  the  observations 
were  obtained.  It  must  be  noted  that,  in  practice,  the  observed  right  aseension 
and  declination  arc  obtained  in  an  astrometric  system,  which  is  defined  by  the 
coordinates  given  to  a  number  of  stars.  Thus  the  difference  between  a  celestial 
and  an  astrometric  system  will  be  due  to  the  errors  in  the  astrometric  system. 
With  this  understanding,  we  will  designate  the  celestial  system  as  the  reference 
system . 

After  transformation  of  the  satellite  and  station  coordinates  into  the  required 
celestial  system  (see  Section  III-F),  the  estimated  right  ascension  and  declina¬ 
tion  can  be  obtained  from  the  following  relationships. 


6  arc  tan  - 

r 

X 


(126) 

(127) 
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where  the  subscript,  0,  indicates  satellite  coordinates  in  the  celestial  system 
and  subscript  sQ  indicates  station  coordinates  in  the  celestial  system. 


(12a) 


The  partial  derivatives  are  obtained  from  the  above  equations  giving 


l  ,  \ 

9x 

-  i 

A  i 

^2 1 

'  \) 

do 

9 

r“ 

£ 

£ 

( x~  -  X  1 

9y 

X 

0 

12 

22 

\  0  so)_ 

do  _  do  do  _  dot 
dz  9x  dy  d  i 


05 

9x 

-  \)  (x»  •  %) 

96 

oy 

"7  WT 

r 

-  \)  (yo  -  \) 

rT 

c 

96 

rv 

9z 

L  xo 

96  _  36  _  96_  _  q 
9x  dy  dz 


where 


[l-]* 


6  0S. 

*11  *12 


13 


6-22  ^23 

•^31  ^32  ^33  J 


T  T 

=  PST  P  N 


(129) 


(130) 


(131) 


(3  32) 


is  the  transformation  matrix  from  the  true  equinox  and  equator  of  date  to  the 
mean  equinox  and  equator  of  the  particular  celestial  system  (see  Section  III-F). 
The  radius 
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■T  ‘[rx  (: 

L  o  \ 


zo  -  z*o 


(133) 


The  matrix  L  normally  will  be  close  to  a  unit  matrix,  and  in  many  cases, 
the  multiplication  of  the  partial  derivatives  by  this  matrix  may  not  be  necessary. 

The  elevation  and  azimuth  rates  can  be  obtained  by  differentiation  of  the  ex¬ 
pressions  for  E  and  A  with  respect  to  time. 


z  r  -  z  r 

T  XT  T  XT 


(134) 


xTyT  -  yTxr 
2  2 
xT+yT 


(135) 


and  the  partial  derivatives  of  the  elevation  rate  with  respect  to  the  orbital  ele¬ 
ments  are 
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(137) 


Similarly,  the  partial  derivatives  of  the  azimuth  rate  with  respect  to  the 
orbital  elements  are 
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where 


XT 


(4^4f 

x^x  +  y  y 
T  T  JT  T 


XT 


XT 


and  the  matrix  S  has  been  defined  previously  (see  Section  III- F) . 

Next  we  will  consider  the  expressions  for  range,  range  rate,  and  range 
acceleration,  which  are,  respectively 


(140) 


rT  =  j-  (xT  x^.  +  yT  yT  +  zT  z^.) 

T 

=  j-  [(X  -  xs)  (x  +UE  ys)  +  (y  -  ys)  (y -u£  xs)  +  (z  -  ze)zj  (141) 

(*2-“e  >f)+(x-xs)(k'-“Ex.)  + 
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(y2  -w2,  Xp  4  (y  -  ys)  (y  -<o?  ys) 


•  2  /  v  •  •  *  2 

+  z“  +  (z  -  zs)  z  -  r^. 


(142) 


Where  w  is  the  rate  of  Earth's  rotation,  wK  -  0.000072921150  radians  per  mean 
solar  second.  The  subscript  s  indicates  station  coordinates. 

The  partial  derivatives  of  range  with  respect  to  the  orbital  elements  are 
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and  the  corresponding  derivatives  of  the  range  rate 


(144) 


- 1 

1  07 

1 

ax 

ar 

T 

X 

W 

rT 

dJr 

dz 

1~ 

- 1 

Qjl  QJ 

] 

af 

T 

JL 

oy 

rT 

af 

T 

az 

(x  +  wE  ys)  -  —  (x  -  Xs) 


(y-w£  xs)  (y-  ys) 


Z - (Z  -  7  ) 

rr 


X  -  xc 


y  -  ys 


z  -  z 


(145) 


(146) 


54 


Finally,  the  derivatives  of  the  range  acceleration  are: 
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VIII.  EQUATIONS  OF  MOTION 
A.  ASPECTS  OF  THE  MATHEMATICAL  MODEL 


The  mathematical  model  of  the  dynamical  system  is  expressed  by  the 
equations  of  motion  with  terms  representing  the  various  forces  acting  o-.  the 
satellite.  The  main  forces  acting  on  a  close  Earth  satellite  are  due  to  Earth's 
gravitational  field,  atmospheric  drag,  Sun's  and  Moon's  gravitational  fields  and 
solar  radiation  pressure.  The  effect  of  these  forces  can  vary  considerably, 
depending  on  the  particular  orbit.  For  example,  a  circular  orbit  of  an  Earth 
satellite  at  800  km  altitude  will  be  perturbed  by  a  maximum  of  about  i.  5  km 
due  to  the  second  harmonic  in  the  Earth's  gravitational  potential,  while  the 
maximum  perturbation  due  to  the  combined  gravitational  attraction  of  Moon  and 
Sun  will  be  less  than  e  meter.  The  perturbations  due  to  atmospheric  drag 
and  solar  radiation  pressure  may  be  of  the  same  order  of  magnitude  at  this 
altitude,  while  at  low  altitudes,  the  effect  of  drag  will  provide  the  principal 
force  .  Therefore,  for  a  close  Earth  satellite,  the  perturbations  due  to  the 
gravitational  attraction  of  the  Moon  and  the  Sun,  and  the  perturbations  due  to 
the  solar  radiation  pressure  may  be  computed  considering  only  first-order 
effects. 

The  equations  of  motion  lor  a  close  satellite  are  written  in  the  previously 
defined  basic  axes  system. 


x  -  +X  +X  +X„  +  X 

G  O  G  C  SR 

(149a) 

y  =Y  +Y  +Y  +Y  +Y„d 

J  G  D  O  I  SR 

(149b) 

z  —  Z  +  z  +z  +  Z  +  Z„ 

G  D  O  z  SR 

(149c) 

The  terms  with  subscripts  G,  D  ,  o  ,  c  and  SR  represent  the  components  of 
acceleration  due  to  Earth's  gravitational  field,  atmospheric  drag,  Sun's 
attraction,  Moon's  attraction  and  solar  radiation  pressure,  respectively.  For 
accurate  computations,  the  independent  time  argument  in  the  equations  of 
motion  must  be  in  a  uniform  time  scale,  such  as  the  atomic  time  (see  Chapter 
H>. 


We  will  now  develop  the  terms  in  a  form  convenient  for  integration  in  rec¬ 
tangular  coordinates. 


B.  EARTH'S  GRAVITATIONAL  FIELD 

The  mathematical  representation  of  the  Earth's  gravitational  field  is  ex¬ 
pressed  by  means  of  the  gravitational  potential  function  which  can  be  written 
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u  =  — 

r 


1  ^  Z  (“)  (SiU  <t>)  {°™  COS  m  A  +  Snir,sin  m  A) 


n=lm=D 


(150) 


Cno=  A  ~  Cn,  Sno  -  0 


where  fi  =  G 

Jn,  Cnm,  and  Snm  are  numerical  coefficients,  RE  the  mean  equatorial  radius  of 
the  Earth,  r  the  distance  of  the  satellite  from  the  center  of  the  Earth,  4>  the 
latitude,  and  P™  the  associated  Legendre  polynomial 


TT\ 


Pn  (x)  =  (l-x“) 


om/2dmPn  (X) 


dx 


m 


(151) 


where  Pn  is  the  Legendre  polynominal.  The  longitude  \  is  to  be  counted  posi¬ 
tive  to  the  east  in  this  application. 


The  harmonics  represented  in  the  gravitational  potential  function  are  called 
spherical  harmonics  (Ref.  18).  If  0  <  m  <  n  they  are  called  tesseral  harmonics 
as  a  special  case  of  the  spherical  harmonics.  If  m  =  0,  they  are  called  zonal 
harmonics,  and  if  m  =  n,  they  are  called  sectorial  harmonics.  The  gravitational 
potential  for  bodies  with  spherical  symmetry  can  be  expressed  by  the  zonal  har¬ 
monics  only,  i.e. ,  the  potential  is  a  function  of  latiiude  and  independent  of  longi¬ 
tude.  For  bodies  of  arbitrary  shape,  the  potential  must  include  the  tesseral 
harmonics,  which  are  dependent  on  both  latitude  and  longitude. 


The  Legendre  polynomials  are  computed  from  the  general  expression: 


1  •  3  •  5  •  •  •  (2n  -  1)  i/z\n  n  (n  -  1)  /z\ 

1  •  2  ■  3  •  •  •  n  )  \  r  /  2  (2n  -  1)  ^r,/ 


n(n-l).(nr2)(n:.3L  ( 

2  •  4  •  (2n  -  1)  (2n  -  3)  \r/  | 


(152) 


The  components  of  acceleration  due  to  the  gravitational  potential  are 


v  9U 
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(153a) 

(153b) 
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(153c) 


Before  we  differentiate  the  potential  function,  we  will  write  U  in  the  following 
form 


u  = 

-  +  U  +  II 
r  z  t 

(154) 

where  n/r  represents  the  potential  of  the  total  mass  concentrated  at  the  center 
of  the  body.  The  potential  functions  representing  the  zonal  and  tesseral  har¬ 
monics  are  Uz  and  UT,  respectively.  By  differentiating  p/r 
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where 
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2  2  2 
=  x  +  y  +  z 

and 


dr  _  x  3r  _  y  8r  _  z 
8x  r  ’  3y  r  *  dz  r 

Now,  computing  the  Legendre  polynomials  and  substituting  into  the  potential 
function, we  obtain  for  the  first  nine  zonal  harmonics 

»,-;44)'  [*(«’  -] 

4  (4)'  [•(?)’  -5] 

4(4)  [“(;)'  -■•!;)'•>] 

44)  *i>f] 
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The  partial  differentials  of  the  above  function  with  respect  to  x,  y,  z  can  be 
written  in  the  form 


(157a) 

(157b) 

(157c) 


where  l designates  the  partial  derivative  of  Uz  with  respect  to  r  when  z 

\  ft  /dU^N 

is  kept  constant,  and  l j  the  derivative  with  respect  to  z  when  r  is  kept  con¬ 
stant.  \  z  h 


By  performing  the  differentiations  we  obtain 


-jlfry  [(-«©'*§]'. 
•fry  -mh 
•fry  t(-  (o'  •  (fy  • 
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and 


V 


Comparing  the  two  derivatives  and  designating 


+  1.5 


f  =-86.625 

5 


f  =  -187.6875 

6 


125- 
r 


f2  =  -7-5(f) 

f3= -17.  5(f)3.  7.  5 1 
f4  =  -39.375  (|)  +  26.  25  (|)2  -  1.875 

(r)5  +  78-75(f)3  '  13' 

+  216.  5625  (^)4  -  59.  0625  +  2. 1875 

f?  =  -402.  1875 (~)7  +  563. 0625  -  216.  5625  19. 

fB  =  -854. 648437  +  1407.  65625  -  703.  828125  (|) 

+  108.  28125  (I)2  -  2.4609375 
fg  =  -1804. 25781  (^)9  +  3418.  59375  -  2111. 48437 

(0 


6875 

A 


+  469.218751 


27.0703125  z 
r 


we  can  write 


r2  L 


M  f2  J2  +  (M  f3  J3  ■  •  •  •  +  i^-)  f9  J9 


(158) 


61 


>1  |N 


. (Sh 


(159) 


and,  designating  the  two  functions  in  the  brackets  as  F(  and  F2  ,  respectively, 
we  can  express  the  accelerations  due  to  the  zonal  har_nonics  as 


X  -  _£*  F 

Xz.  "  3  F] 

r 


— —  -  y  a  -HZ  F 
ay  *  -3  'i 


^  =  z,  =  -7(f  F,  -0 


(100a) 


(160b) 


(160c) 


We  will  now  develop  the  tesseral  harmonics  considering  the  gravitational 
potential  up  to  the  fourth  order  harmonics. 


U 

t  r 


r A  /R  \n  m 

^  VT/  <Sln  ^  <Cnm  C°S  m  X  +  Snin  Sin  m  X)  <161) 

n  =  2  rn=l 


For  integration  in  rectangular  coordinates,  it  is  convenient  to  express  the 
trigonometric  functions  in  terms  of  x",  y",  z",  which  are  the  satellite  coordinates 
in  an  Earth  fixed  terrestrial  system.  The  integration  itself  is  done  in  the  basic 
system  defined  in  Section  III-B.  The  required  coordinate  transformations  are 
given  in  Section  III.  To  facilitate  writing  we  will  adapt  the  notation 

x  =  x" 


y  =  y” 

z  =  z” 

With  this  notation,  the  trigonometric  functions  can  be  expressed  as 
sin  4>  =  ~ 
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where 


r2  =  x2  +  y2  +  zJ 

The  associated  Legendre  polynomials  then  are 


2  +  V2 

P  =3 

2  r2 


1  _  1  (5  z2  -  r2)(x2 


1/2 


P3  =2 


=  15liLiXl 

r3 


P_  =  15 


/  2  2\3/2 

(x  +y  . 


p;  = 


=  5  z  (7  z2  -  3  r2)  (x2  +  y2) 
2  4 


1/2 


p2  -  15  (7  z2  -  r2J_(x2  ->■  y2j 
4  2  4 


-  (x2  udl 


3/2 


P  =  105 

4  4 

r 


P  =  105 


1*1 


Similarly,  we  can  express  the  longitude  dependent  functions  in  terms  of  the 
rectangular  coordinates. 

Then  we  can  write 


U  =  n(f  +  f  +  f  +  f  +  f  +  f  +  f  +fj 

T  '  22  31  32  33  41  42  43  44 


(162) 
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wh'4; 


'R 


U2  =\*fj  ~  [C,2  (X2  -  y2)  +  S22  (2xy)] 


si 


■i  ■> 


K  j  3  (r>z“  -  r~) 
2r‘ 


[c„  *  *  s3,  y] 


=fi)  [c32  <*2  -  y2> +  S-J2  M 


f„  =(t)  T  [°33  X  <x2  ■  3y2)  *  S33  y  ,3x2  •  y2)] 
\  /  r 


•■41 


\)  5z  (7z2  -  3r2) 

c 

2r 


[c4i  x  +  s4i  y] 


f  =fV\  15  (7z2  -  r2) 
42  "  V  r  /  2r5 


R 


[c42  (X2  -yVS42  (2xy)] 


f33  =(-#)  ~  [C43  X  (X"’  -  3y2>  *  S43  y  <3X  -  y2>] 


=  vf)  ^  [C44  ,X‘  ‘  6x2  y2  -  yJ  >  ‘  S44  4xy  <x2  '  y2  >] 


Performing  the  differentiations  we  obtain 
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-f~{C22  X  l2*  ~  5  "  y2>]  +  S22  2y  (I*2  "  5X2>} 


3  m-  R 
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15  m  z 
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{c32  x  [2r2  -  7  (x2  -  y2)]  +  S  ,  2y  (r2  -  7x2)} 
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15nR„ 


~{c.i3  [  3r2  (x2  -  y2)  -  7x2  (x2  -  3y2)]+  xy  [er2  -  7(3x2  -  y“ )]  J 


5*i  R  .  z 


2r 
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+  S41  21xy  (r2  -  3z2)} 
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315*i  R  z 

E 
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S44  4y  [r2  (3x2  -  y2)  -  9x2  (x2  -  y2)]} 


(163) 
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- (3r“-7z2)  [C31x+S31y] 
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315k  R  z 


(r2  -  3z2  )  [C42(x2-y2)  +  S42(2xy)] 


105k  R 


E  (r2  -  9z2 )  [C43  x  (x2  -  3y2 )  +  S43  y  (3x2  -  y2)] 


945k  R  z 

E 


[c44  (X4  -  6x2y2  +  y4)  +  S4>,  4xy  (x2  -  y2)] 


(165) 


The  accelerations  due  to  the  gravitational  potential  then  can  be  written 


X  _  _  M-  x  ^  x  x 


(166a) 
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til  ,  Y  ►  Y 

r3  /■  T 

(166b) 

—  +  z  +  z 

3  Z  T 

(166c) 

C.  ATMOSPHERIC  DRAG 

The  acceleration  due  to  drag  on  a  satellite  is  a  function  of  atmospheric 
density,  p,  relative  velocity  of  the  satellite  with  respect  to  the  atmosphere, 
vr,  satellite  mass,  m,  drag  coefficient,  CD,  and  reference  area,  A. 

The  atmospheric  density  can  vary  considerably  and  is  rather  difficult  to 
evaluate  with  an  accuracy  necessary  for  precise  orbit  determination,  particularly 
at  lower  altitudes  where  its  effect  is  large.  It  is  a  function  of  altitude,  exospheric 
temperature  (e.  g.  Ref.  19),  and  the  relative  position  of  the  Sun. 


If  the  density  is  given,  the  expression  for  the  acceleration  is 
pv2  C  A 

_  K  r  p 

ao  2  m 


(167) 


The  relative  velocity,  vr,  can  be  obtained  on  the  assumption  that  the  atmos¬ 
phere  rotates  with  the  same  angular  velocity  as  the  Earth.  The  relative  velocity 
components  in  the  basic  axes  system  then  are 


v  =  x  +  OJ  V 
vx  A 

(168a) 

vy  =  y- ^Ex 

(168b) 

V  =  z 

z 

(168c) 

Where  is  the  angular  velocity  of  the  Earth's  rotation. 

Then 

222 

V  =  (V2  +  V2  +  V2  ) 

(169) 

and  the  components  of  the  accelerations  in  the  basic  axes  system  are 


XD=-Pvr 


v*  2m 


(170a) 
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Pv  V 

r  v 

C  A 

[> 

( L  70b) 

2m 

PV,  V, 

C»A 

(170c) 

2  m 

D.  SUN’S  GRAVITATIONAL  ATTRACTION 

The  perturbing  action  of  the  Sun  on  a  eiose  satellite  orbit  can  be  expressed 
by  tiie  so  called  disturbing  function.  II  the  ratio  of  the  satellite  distance  from 
the  center  of  Earth  to  the  Sun's  distance  is  sufficiently  small,  the  disturbing 
function  can  be  expanded  in  power  series  whieh  converge  rapidly.  The  disturbing 
function  is  given  in  a  general  form  as 


xx  +  yy  +■  zz 

J” 


(171) 


where 


M-  ,=  GM  , 


and  x,  y,  z  and  x  ,  y  ,  z  are  the  satellite  and  Sun's  coordinates,  respectively. 
Ais  the  distance  between  the  satellite  and  the  Sun, 


A"  -  (x  -  x  )-  My  -  y  >)2  +  (z.~  zf 


(172) 


The  listurbing  function  can  be  expanded  in  powers  of  r/r,  (Ref.  20)  giving 


P; 

r 


(f*)  (l  eosi0  '  2)  +  (r^  (l  eos2  0  "  |  cos  °)  4  •  •  •  (173) 


where  the  functions  in  the  brackets  are  the  Legendre  polynomials  which  were 
treated  in  Chapter  VIII,  Section  B. 


0  is  the  angle  between  the  directions  to  the  satellite  and  the  Sun. 
Consequently, 

xx  +yy  +zz 

cosO  — -r - -  <m) 

rr 
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and  we  can  write  the  disturbing  function  neglecting  all  higher  order  terms  as 

R  =~3  ~~  (xx0  +  yy0  +  zzo)2 "  i (x2  + y2  +  (175> 

ro  L2ro 

The  accelerations  due  to  the  Sun's  attraction  then  are  obtained  by  differentiation. 


8R 
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3r  f  nr?0  +  ??o)-  x 
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(176a) 
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ay  " 

_rr 

3r  ria  (?fQ+  r?i7Q+  y 

(176b) 
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[3r  ?0  (??0+  mo+  zj 

(176c, 

where  r? f  £,  and  £  ,  rjy  are  the  direction  cosines  of  the  satellite  and  the 
Sun's  radius  vector,  respectively. 

The  direction  cosines  of  the  Sun  can  be  obtained  with  sufficient  accuracy 
from  the  following  equations. 


?  =  cos  vc. 

(177a) 

i73  =  sin  vQ  cos  € 

(177b) 

-  sin  wsin  e 

(177c) 

where  the  Sun's  longitude,  v0  ,  from  the  mean  equinox  of  date  is 
v3  =  4. 8883394  +  0.  017202791  d 

+  0.  03345  sin  (6. 2482906  +  0.  0172019697  d)  (178) 

and  the  mean  obliquity  of  the  ecliptic,  e  ,  is  obtained  from 

€  =  0.  40920619  -  0.  6218433  x  10**  d  (179) 

The  Sun's  distance  from  the  center  of  Earth  can  be  obtained  in  kilometers  from 
r0  =  1.  496  x  108  [l  -  0.  016725  cos  (6.  2482906  +  0.  0172019697 
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These  expressions  are  based  on  the  fundamental  ephemerides  as  given  in  Ref.  1. 
The  angles  arc  in  radians  and  d  is  in  ephemeris  days  from  1950  January  1,  0h  E.T. 


E.  MOON’S  GRAVITATIONAL  ATTRACTION 


Since  the  Moon  ;s  much  closer  to  the  Earth  than  the  Sun,  the  convergence 
of  the  series  in  the  disturbing  function  will  be  slower.  Therefore,  for  higher 
altitude  satellites  and  accurate  orbit  computations  it  may  be  necessary  to  con¬ 
sider  more  terms.  To  a  first  order  accuracy,  we  can  write  the  accelerations 
as  obtained  in  the  previous  section. 
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3r  (??c+  11 n:  +  x 
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_3r  (?et+  m*  SfJ-  * 


(181a) 


(181b) 


(181c) 


Where  ?,  17,  %  arc  defined  in  the  previous  section,  and  are  the 

direction  cosines  of  the  Moon's  radius  vector  from  the  center  of  the  Earth. 
They  can  be  computed  from  the  following  expressions 


-  cos  A  cos  Q  -  0. 99596  sin  AT  sin  Q  (182a) 

Tl.=  cos  e  (cos  At  sin  Q  +  0. 99596  sin  ATcos  ft) 

-  0. 08976  sin  €  sin  AT  (182b) 

?,.=  sin  e  (cos  AT  sin  Q  +  0.  99596  sin  \T  cos  ft) 

+  0.  08976  cos  e  sin  XT  (182c) 

where  e  is  the  mean  obliquity  of  the  ecliptic  given  by  Eq.  (179). 


^is  the  angular  position  of  the  Moon  measured  in  the  Moon’s  orbital 

plane  from  the  ascending  node  on  the  ecliptic,  and  is  obtained  from  the  ex¬ 
pression 
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0.  91215551  +  0.23089572  d  +  0. 10980098  sin  <3.  7617316 

+  0.22802714  d)  (183) 

S7  is  the  longitude  of  the  ascending  node  of  the  lunar  orbit  on  the  ecliptic, 
measured  from  the  mean  equinox  of  date. 

ft  =  0.21140807  -  0.  0009242193  d  (184) 

The  lunar  distance  can  be  computed  from 

r,  =  384400  |l  -  0.  054900489  cos 

In  the  above  equations,  the  angles  are  given  in  radians,  r,  is  in  kilometers, 
and  d  is  defined  in  Section  VIII-D.  The  higher  order  terms  have  been  omitted 
but  the  equations  are  sufficiently  accurate  for  our  purposes. 


(3.  7617316  +  0.  22802714  d)  (185) 


F.  SOLAS,  RADIATION  PRESSURE 


The  energy  of  solar  radiation  which  is  imparted  to  a  celestial  body  per  one 
square  meter  of  surface  area  per  one  second  of  time  is  (Ref.  21) 


S  =  1350 


ft); 


joules 


(186) 


m  sec 


Where  r0  is  the  mean  distance  of  Earth  from  Sun,  and  r  is  the  distance  of  the 
body  from  Sun. 


The  force  on  the  body  generated  by  this  energy  can  be  expressed  as 

F  ,SJl_tR)  2a  ["newtons'] 
sr  c  L  J 


(187) 


where  R  is  the  reflection  coefficient,  R  =  0  for  an  absolute  black  body,  and 
R  =  1  for  a  specular  surface,  c  is  the  velocity  of  light,  and  a  is  the  angle  of 
incidence  with  the  surface. 


For  a  sphere,  this  equation  reduces  to 

newtons 


_2S(H  R) 
SR  C 


m 


(188) 


The  acceleration  produced  by  this  force  will  be 


FI  mm 

SR  SR  m 


m 

2 

sec 


(189) 
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where  A  is  the  reference  area,  and  m  the  mass  of  the  body. 

The  direction  of  the  acceleration  is  away  from  the  Sun,  and  thus  it  opposes 
the  Sun's  gravitational  attraction.  The  components  of  the  acceleration  are 
found  by  multiplying  aS[<  by  the  direction  cosines  £  ,  n  .  C 


XSH 

*  aSR  ^  ' 

(190a) 

Ys,« 

=  *  aSK  rl 

(190b) 

ZS,, 

=  -  a.  ? 

SR  ' 

(190c) 

It  must  be  noted  that  the  direction  cosines  ?Q,  r)Q,  ^  are  defining  the 

direction  Earth- Sun.  However,  for  close  satellites,  the  error  introduced  by 
using  them  for  satellite-Sun  direction  is  negligible. 

The  radiation  pressure  is  experienced  by  the  satellite  only  at  times  when 
it  is  in  the  sunlight  The  angle  0  between  the  two  vectors  in  the  directions  to 
satellite  and  Sun  is  found  from 

cos  0  -  J-S  m  +  cf  (191) 

The  satellite  is  always  in  sunlight 
when 


r  sin 0  >  R,. 

r. 

If  r  sine  <  Rt.,  the  satellite  is  in  sunlight  only  when  simultaneously  cos  0  >  0. 

G.  INTEGRATION  OF  THE  EQUATIONS  OF  MOTION 

Integration  methods  can  be,  generally,  placed  in  three  groups:  (1)  self¬ 
starting  techniques ,  (2)  difference  techniques,  and  (3)  predictoircorrector  tech¬ 
niques.  Although  all  three  techniques  have  been  used  in  orbit  integration,  the 
self- starting  methods  have  a  distinct  advantage  in  orbit  determination  programs. 
The  advantage  stems  from  the  fact  that  the  integration  must  be  done  between 
two  arbitrary  observation  times,  but  the  last  two  techniques,  normally,  employ 
a  constant  integration  interval.  Moreover,  the  differential  equations  for  a  close 
Earth  satellite  are  of  second  order  in  which  the  first  derivatives  are  included. 
Most  of  the  standard  methods  of  the  second  and  third  type  used  in  astronomical 
orbit  integration  have  been  developed  for  the  special  case  in  which  the  first 
derivatives  are  absent  (Ref  22).  The  disadvantage  of  the  self- starting  methods 
lies,  mainly,  in  the  computing  time  required  for  the  integration.  However, 
this  can  be  compensated  to  a  large  extent  by  incorporating  methods  which  vary 
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the  integration  step  size  and/or  estimaU  and  adjust  the  error  according  to  some 
criterion. 


Errors  in  numerical  integration  are  due  to  two  causes:  truncation  and  round¬ 
off.  Truncation  errors  are  caused  by  the  replacement  of  the  actual  differential 
equations  by  ordinary  difference  equations  and  omission  of  higner  order  terms. 
The  round-off  errors  are  due  to  the  finite  number  of  digits  carried  in  the  com¬ 
putations.  Not  much  can  be  done  about  the  round-off  errors  besides  increasing 
the  number  of  digits.  The  truncation  errors,  however,  can  be  controlled  to 
a  certain  extent  by  varying  the  integration  interval.  Two  types  of  techiques 
generally  are  used  for  this  purpose.  One  technique  estimates  the  truncation 
error  and  adjusts  the  integration  interval  and  error  after  a?i  initial  integration 
step  has  been  computed.  This  involves  additional  computations  and  thus,  in 
part,  cancels  the  benefits. 

In  cases  where  the  nature  of  the  differential  equations  is  known,  the  criterion 
can  be  determined  in  advance  and  each  step  size  determined  automatically. 

Thus,  considerable  savings  in  time  can  be  gained.  A  technique  of  this  kind  is 
used  in  the  present  application.  The  rule  for  the  step  size  is 

Ah  -•  k  r  (192) 

where  Ah  is  the  integration  interval,  r  is  the  distance  from  the  center  of  the 
Earth,  and  k  is  a  constant  for  the  orbit,  determined  as  follows. 


At  perigee 


and 


thus 


r  A 9  =Vper  Ah 


Ah  =~~  r  =  k  r 

Vper 


A  9-  /CY'2 

'  1  +  e 


(193) 


Where  A9  is  dependent  on  the  integration  method.  It  was  found  that  optimum 
values  of  A9  are  obtained  from  the  following  empirical  expression: 


A9  =  Cj  +  c2  exp 


(194) 


Where  e  is  eccentricity  and  c  x ,  c2,  c3  are  constants  which  depend  on  the  in¬ 
tegration  method.  The  criterion  gives  approximately  equal  arc  lengths  for  the 
integration  intervals. 
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The  most  popular  of  the  self- starting  techniques  has  been  the  fourth  order 
Runge-Kutta  method.  A  similar  method,  which  compares  rather  favorably 
with  the  Runge-Kutta  method,  has  been  developed  by  Bowie  (Ref.  23).  This 
method  is  based  on  the  assumption  that  the  second  derivatives  vary  quadratieally 
over  the  integration  interval. 

In  the  present  effort,  these  methods  are  applied  in  a  direct  integration  mode. 
In  a  broader  sense  it  could  be  classified  as  the  Cowell's  method.  A  different 
principle  is  employed  in  the  Enekc's  method.  In  this  method  the  coordinates 
are  not  obtained  directly,  but  rather  the  integration  is  performed  on  the  differ¬ 
ence  between  the  actual  coordinates  and  the  coordinates  of  a  conie  section. 

The  conic  section  is  obtained  from  the  position  and  velocity  components  at  a 
particular  instant,  called  the  epoch  of  osculation.  The  departures  from  the 
osculating  orbit  are  called  perturbations.  If  the  perturbations  are  small,  they 
can  be  expressed  by  fewer  significant  figures,  which  permits  larger  intervals 
than  with  the  direct  integration  methods.  When  the  perturbations  increase 
to  an  intolerable  size,  a  rectification  of  the  orbit  is  required.  The  position 
and  velocity  components  are  determined  at  a  new  epoch  and  the  integration 
restarted. 

The  application  of  the  Eneke's  method  in  close  satellite  orbit  determination 
programs  is  impaired  because  of  the  irregular  intervals  between  observations 
which  complicates  the  solution  of  the  two-body  equations.  In  any  case,  a  definite 
superiority  between  the  Eneke  and  Cowell  type  methods  in  modern  computer 
applications  has  not  been  established. 
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APPENDIX 


A.  EQUATIONS  OF  MOTION  IN  RECTANGULAR  COORDINATES 

Equations  will  be  developed  for  motion  in  an  inverse -square  central  force 
field.  According  to  Kepler’s  laws,  the  path  of  such  motion  is  a  conic  section 
with  the  central  mass  occupying  one  focus.  We  will  first  consider  an  elliptic 
motion.  An  eliipse  is  formed  if  a  point  moves  in  such  a  manner  that  the  sum  2a 
of  its  distances  r x  and  r2  from  two  fixed  points,  the  focuses,  is  a  constant,  i.e. , 

i’j  -  r,  =  2a  =  const  (A.l) 

where  a  is  the  semimajor  axis.  In  a  Cartesian  axes  system  with  the  origin  at 
one  focus 

x2  +  y2  =  r2  =  r2  (A .  2) 

and 
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where  b  is  the  semiminor  axis.  Substituting  (A.  6)  into  (A.  5)  we  obtain 


■> 

“  -  C  (A.  7) 

which  is  an  expression  for  the  semilatus  rectum.  Equation  (A. 4)  can  now  be 
w  r  it  ten 

c  cos  u5  x  +  e  sin  3  y  *  C  =  r 

where  c  is  the  numerical  eccentricity  (or  simply  eccentricity) 


c  = 


£ 

a 


Substituting  A  =  e  cos  3,  and  B  =  c  sin  3  we  obtain 

Ax  +  By  +  C  ^  r  (A.  8) 

which  is  the  equation  of  an  ellipse  in  rectangular  coordinates.  Differentiating 
Eq  (A.  8)  with  respect  to  time  we  obtain 

(A  r  -  x)  x  +  (B  r  -  y)  y  =  0  (A .  9) 

For  a  motion  in  a  central  force  field  with  the  force  varying  inversely  as  the 
square  of  the  distance  from  the  central  mass,  the  following  relationship  holds 

h2  -  m  C 

where  h  is  the  angular  momentum  about  the  center  of  force.  Expressed  in  rec¬ 
tangular  coordinates,  the  angular  momentum  is 

h  =  x  y  -  y  x 

which  gives 

(x  y  -  y  x)2  =  P-C  (A.  10) 

Equations  (A. 8),  (A. 9)  and  (A.  1C)  completely  define  the  motion.  Identical  equa¬ 
tions  can  be  obtained  for  a  parabola  and  a  hyperbola  following  the  same  procedure. 


From  the  definition  of  A  and  B  it  follows  that 
e  =  (a"  +  B2)1/2 

For  a  case  y  =  0,  the  expressions  reduce  to 
(1  -  A)  x  =  C 
(1  *  A)  xu  =  B 

x.  t  (»*  c >1/2 


(A.  11) 


(A.  12) 
(A.  13) 
(A.  14) 
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Solving  for  the  planar  orbital  elements  A,  B,  and  C 


A  = 

B  - 

C  - 


x  x  y 

w.  W  **  U- 

M- 

fT~ 


(A.  15) 
(A.  1(5) 
(A.  17) 


Equations  (A.  8),  (A.  9)  and  (A.  10)  can  be  solved  for  the  velocity  components 


(A. 18a) 

(A. 18b) 


where  the  upper  signs  are  for  a  direct  motion,  and  the  lower  signs  are  for  a 
retrograde  motion.  When  y  =  0,  x  =  x^  =  r  and  the  equations  reduce  to 


x 


(A.  19) 


'(A  -  1) 


(A.  20) 


B.  EQUATION  OF  TIME 

We  can  write  Eq  (A.  8)  in  a  canonical  system  in  which  the  x-axis  points  to¬ 
ward  the  pericenter  and  the  coordinates  are  normalized  dividing  by  C.  In  this 
system  A  =  -  e ,  B  =  0 ,  and  C  =  1 . 

/  2  2V 

-c  x  +  1  =  ^x‘+  y  )  (A. 21) 

This  equation  can  be  solved  for  y 
r  2  _2  _  1 

y  =  ±  [(e  -  1)  x  -  2e  x  +  ij  (A.  22) 

or  for  x 


x 


(A.  23) 


where  only  a  negative  sign  before  the  radical  applies  in  the  chosen  axes  system. 
The  relation  between  the  coordinates  x,  y  and  x  ,  y  is  given  by 


C  x  -  x_  *  y_  ria 

Cy  -  x_Pa  -  y_  §a 


(A.  24a) 
(A.  24b) 


The  normalized  area  swept  out  by  a  radius  vector  from  the  x-axis  to  a  point 

P  (x  .  y  )  is  obtained  by  integration 

k  k  k 

V  \  —  — 

r  p  y. 


I  k.  * 


and  substituting  y 


s,  )  [<«* 


1)  x2  -  2e  x  +  1 


.i/2 

.J  d  x 


x,  y, 

k  k 


(A.  25) 


Considering  that 


kP  ~  1  +  e 


The  integral  (A. 25)  can  be  obtained  for  three  cases 


2  3/2 

2(1  -  e“) 


j|  -  arc  sin  [e  +  (1  -  e2)  xk]  -  e(l  -  e2)17"  yk  |  if  e 


(A.  26) 


t  3/2 

Z(e“  -  1) 


if  e  >  1 


| e(e2  -  1)1/2  yk  +  In  [e  -  (e2  -  l)l/2  yk  -  (e2  -  1)  xj| 


(A.  27) 


y,3  x  y 

*  «  k  -e  , 

^*32  >  if  6  -  1 


(A.  28) 


From  Eq  (A.  10)  we  can  obtain  the  normalized,  constant  area  rate 


n = 


(A.  29) 
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Returning  to  Eq  (A.  26),  it  can  be  shown,  utilizing  the  geometry  in  Fig.  A-2,  that 
sin  0k  =  e  +  (x  -  e2)  xk 
and 

2  1/2  - 

cosek  =  (1  -  e  )  yk 

Thus  0k  is  the  complement  of  the  eccentric  anomaly  Ek 


Ek  ,|-ek  if  y  2.0 

Ek  *  I1  ‘ek  ify<0 

and  we  can  write 

2  — 

cos  E  =  e  +  (1  -  e  )  x 

K  K 

o  1/2 

sin  Ek  =  (1  -  e“)  yk 

which  shows  that  Eqs  (A.  26)  and  (A.  30)  are  another  form  of  Kepler's  equation 
Mk  =  Ek  -  e  sin  Ek  (A. 31) 

where 


Kepler's  equation  is  transcendental  in  E  and  can  be  solved  by  successive  approxi¬ 
mations,  utilizing  the  derivative 


dtk 


_ 1 _ 

Kc(l  -  e  cos  E  ) 


(A.  32) 


Similar  considerations  apply  to  the  hyperbolic  case.  Differentiating  Eq  (A. 27) 
we  obtain 
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I  9  ’/2 

Bk  e  -  (Q~  -  l)  Xk-  (e2  -  1)  “  yk 

For  a  parabola,  a  closed  form  solution  can  be  obtained.  Given  time  t^  from 
pericenter  to  P  (x, ,  y  ),  the  corresponding  normalized  area  is 

K  K 

A  -  (A.  34) 

*  K  K. 

C 

where 


If  e  =  1,  Eq  (A,  21)  gives 


Substituting  Eq  (A.  35)  into  (A.  28)  we  obtain  a  cubic  in  yk 


(A. 35) 


V3*3y  -  12A  =  0 

Jk  k 

The  discriminant  of  this  equation,  A>  0,  and,  therefore,  it  will  always  have  one 
real  root 

yR  =  [6Ak  -  (36Ak  +  1)U2 ]  +  [6Ak  -  (36Ak  +  1)1/2J  (A,  36) 

and  xk  is  obtained  from  Eq  (A. 35). 
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